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SECTION  I 
INTRODUCTION 


The  Specter  computer  codes,  as  described  in  references  1-3,  determine  the 
radiation  environment  resulting  from  a high-altitude  nuclear  explosion  and  the 
irradiation  of  satellites  in  that  environment.  This  report  describes  recent 
improvements  that  have  been  made  in  the  computation  of  the  distributions  of  the 
nuclear  debris  and  trapped  fission-decay  electrons  in  the  magnetosphere,  in  the 
evolution  of  these  distributions  with  time,  and  in  the  ultimate  trapped-electron 
environment  that  may  ensue  in  the  event  of  a nuclear  exchange. 


A new  electron  injection  model  is  described  in  Section  II.  It  incorporates  a 
model  for  the  fission-debris  distribution  that  is  physically  more  realistic 
than  the  "confined"  and  "jetting"  debris  models  that  were  used  in  the  previous 
versions  of  the  Specter  codes.  The  magnetic  tube  containing  the  debris  now 
convects  outward,  toward  higher  L values,  owing  to  the  enhanced  plasma  pressure 
in  the  tube.  The  trapped  electrons,  injected  by  the  decay  of  the  debris,  drift 
outward  as  well  as  eastward  while  they  reside  in  the  debris  tube.  The 
redistribution  of  the  electrons  associated  with  their  displacement  in  L is 
computed  assuming  conservation  of  the  first  two  adiabatic  invariants  of  the 
electron  motion.  Since  the  lower-energy  electrons  have  lower  eastward-drift 
velocities,  they  spend  more  time  in  the  debris  tube,  and  are  therefore  trans- 
ported to  higher  L values.  Accordingly,  the  energy  spectra  of  the  trapped 
electrons  rapidly  become  softer  toward  higher  L values.  This  result  of  the 
new  injection  model  is  in  agreement  with  the  nuclear-tests  data;  it  constitutes 
a significant  improvement  in  the  computation  of  the  trapped-electron  distribu- 
tion. A comparison  of  the  computed  fluxes  with  the  nuclear-tests  data  is 
described  in  Section  III. 


As  the  population  of  trapped  electrons  increases  due  to  injections  by  multiple 
nuclear  explosions,  the  decay  rate  is  expected  to  increase  due  to  enhanced 

i 

pitch-angle  diffusion.  This  effect  is  taken  into  account  by  incorporating  into  j 
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Specter  an  empirical  model  for  the  time-dependent  loss  rates  based  on  the 
work  of  M.  Schulz  (Rtf.  4).  A discussion  of  this  model  together  with  the 
self-consistent  flux  levels  for  the  saturation  environment  is  given  in 
Section  IV.  Section  V treats  the  conditions  under  which  two  or  more  bursts 
interact  such  as  to  alter  the  injection  efficiency  of  the  individual  bursts. 

The  trapped  electron  loss  model  in  Specter  is  highly  empirical,  and  makes  no 
distinction  between  the  various  processes  that  remove  electrons,  such  as 
pitch-angle  diffusion  at  high  altitudes  and  scattering  in  the  atmosphere. 

There  has  been  no  attempt  made  until  recently  to  include  processes  that 
redistribute  electrons  between  adjoining  L-shells.  The  new  injection  model, 
however,  gives  much  improved  radial  distributions;  accurate  enough  to  warrant 
an  investigation  of  subsequent  processes  that  alter  the  distribution.  Effects 
similar  to  those  that  lift  the  debris  and  electrons  to  high  L-shells  can  also 
result  in  radial  diffusion  of  trapped  electrons.  Radial  diffusion  is  very 
important  for  energetic  particles  in  the  outer  radiation  belts.  Section  VI 
addresses  the  effects  of  radial  diffusion  on  artificially  injected  electrons, 
and  the  corrections  to  the  distribution  that  may  be  necessary. 

The  electron  loss  model  is  always  subject  to  improvement  as  new  data  become 
available.  The  model  is  weak  in  some  regions  where  adequate  data  have  never 
been  taken.  This  is  especially  true  at  low  altitudes,  particularly  in  the 
altitude  region  between  the  local  and  global  cutoffs.  New  data  are  gradually 
appearing,  however,  and  may  soon  warrant  substantial  revisions  to  the  loss 
model.  Section  VII  is  an  examination  of  the  implications  of  new  electron  loss 
data  to  the  Specter  loss  model.  We  review  (in  VII)  the  recent  data  bodies  and 
describe  (In  VII. 2)  the  comparison  of  the  new  injection  model  with  both  old 
and  new  data.  A concluding  subsection  (VII. 3)  addresses  the  anticipated 
improvements  and  alterations  that  are  needed  to  bring  the  model  into  agreement 
with  the  data.  Recommendations  for  further  improvements  of  the  Specter  codes 
are  made  in  Section  VIII. 
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SECTION  II 


NEW  ELECTRON  INJECTION  MODEL 
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I . INTRODUCTION 

In  order  to  compute  the  initial  distribution  of  the  fission  beta-decay 
electrons  which  are  injected  into  the  radiation  belt  due  to  a high-altitude 
nuclear  detonation,  it  is  necessary  to  know  the  spatial  and  temporal  distri- 
bution of  the  fission  fragments.  A new  model  has  been  developed  for  the 
computation  of  the  fission-debris  distribution.  This  model  is  physically 
more  realistic  than  the  "confined"  and  "jetting"  debris  models  that  were 
previously  used  in  the  Specter  Codes. 

Briefly,  the  phenomenology  of  the  ionized  debris  issuing  from  a nuclear 
explosion  above  the  sensible  atmosphere  of  the  earth  is  as  follows.  Initially, 
the  ionized  debris  expands  radially  against  the  geomagnetic  field,  forming  a 
magnetic  bubble  by  excluding  the  field  from  the  highly  conductive  plasma  and 
compressing  it  into  a spherical  shell.  During  the  expansion  and  subsequent 
collapse  of  the  bubble,  the  ionized  debris  and  energetic  air  ions  escape  from 
the  local  entrapment  and  move  toward  higher  and  lower  altitudes  along  the 
earth's  magnetic  field.  The  bundle  of  field  lines  which  contain  the  debris 
and  hot  ionized  air  are  referred  to  as  the  debris  tube.  Through  the  beta- 
decay  process,  which  is  very  rapid  at  early  times,  high  fluxes  of  relativistic 
electrons  are  emitted  by  the  fission  fragments.  A large  fraction  of  the 
electrons  become  trapped  in  the  debris  tube.  Owing  to  the  configuration  of 
the  geomagnetic  field,  the  trapped  electrons  drift  toward  the  east  and  the 
ionized  debris  and  air  drift  toward  the  west.  The  oppositely-drif itng 
particles  produce  a charge  separation  across  the  tube.  As  depicted  in  Figure 
1,  a net  positive  charge  forms  at  the  western  surface  of  the  tube  and  a net 
negative  charge  forms  at  the  eastern  surface.  These  charges  produce  an 
electric  field,  E,  in  the  debris  tube  directed  toward  the  east,  which  causes 
the  debris  to  drift  outward  toward  higher  L values  with  the  I X B drift  velocity. 
Hence,  as  the  debris  drifts  outward,  the  emitted  electrons  drift  eastward 
forming  the  radiation  belt.  The  process  continues  until  the  debris  either  enters 


11 


A 

* 


Figure  1.  Illustration  of  electric  field  across  debris  tube  and  current  tending  to 
neutralize  the  excess  charges. 
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the  atmosphere  because  of  its  motion  along  the  magnetic  field  or  drifts  outward 

beyond  the  boundary  of  the  trapping  region.  The  energy  spectra  of  the  trapped 

electrons  become  softer  at  the  higher  L values  principally  because  the  electrons 

become  displaced  upward  in  L before  they  drift  out  of  the  convecting  debris  tube,  * 

and  the  lower  energy  electrons  reach  higher  L values  because  their  eastward  drift 

velocities  are  smaller. 

The  outward  motion  (E  X B drift)  of  the  debris  depends  sensitively  on  the 
coupling  of  the  magnetosphere  to  the  ionosphere.  The  charge  separation  across 
the  tube  tends  to  be  neutralized  by  currents  in  the  thermal  plasma  which,  as 
shown  in  Figure  1,  flow  along  magnetic  field  lines  at  high  altitudes  and  across 
magnetic  field  lines  in  the  ionosphere  where  the  collision  frequency  is 
sufficiently  high.  The  ionizing  radiations  from  the  nuclear  device  greatly 
increase  the  conductivity  of  the  local  ionosphere.  This  effect  tends  to  dis- 
charge the  tube,  hence  to  retard  the  outward  motion  of  the  debris.  On  the  other 
hand  the  neutralization  current  along  the  magnetic  field  may  become  so  high 
that  "anomalous  resistivity"  may  develop  along  the  field  lines  in  the  upper 
ionosphere.  Such  a resistivity,  which  is  due  to  a plasma  instability  (Refs. 

5-9)  tends  to  decouple  the  tube  from  the  ionosphere,  hence  to  attenuate  the 
neutralization  process. 

In  the  new  injection  model,  the  computation  of  the  outward  motion  of  the 
debris  is  based  on  the  assumption  that  magnetic  field  lines  are  equipotentlals. 

This  assumption  greatly  simplifies  the  problem  because  the  debris  remains  in  a 
magnetic  flux  tube.  The  neutralizing  effect  of  the  ionospheric  current  is 
determined  by  using  an  equivalent  electrical  circuit  to  computer  the  net  charge 
across  the  tube.  This  circuit  consists  of  the  effective  capacitance  and 
inductance  of  the  magnetic  tube  and  the  variable  resistance  of  the  current  path 
through  the  ionosphere.  It  is  driven  by  the  azimuthal  drift  current  due  to  the 
trapped  electrons  and  the  ionized  debris  and  hot  air. 

The  injection  model  takes  as  input  the  velocity  and  spatial  distribution 
of  the  debris  and  hot  ionized  air  escaping  from  the  magnetic  bubble.  These  data 
were  supplied  by  the  Naval  Research  Laboratory.  The  model  then  computes  the 
distribution  of  the  debris  in  the  magnetosphere  by  combining  the  motion  of  the 
debris  along  the  magnetic  tube  with  the  cross  L motion  of  the  tube.  Finally, 
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the  debris  distribution  — specifically,  the  number  density  of  the  fission  fragments 

as  a function  of  B,  L,  and  time  — are  used  to  compute  the  trapped  electron  distri- 
bution, including  the  redistribution  of  the  electrons  resulting  from  the  motion 
of  the  debris  tube. 

2.  MODEL  OF  CROSS-L  MOTION  OF  DEBRIS 

a.  Equipotential  Field  Lines 

Owing  to  the  highly  complex  behavior  of  a dense,  high-energy  plasma  injected 
into  the  magnetosphere,  some  simplifying  assumptions  must  be  made  to  estimate 
the  motion  of  the  nuclear  debris.  One  approximation  that  is  made  is  that  the 
j magnetic  field  lines  containing  the  debris  are  equipotentials . This  approxi- 

mation is  commonly  made  in  magnetospheric  physics  because  mobilities  of  the 
thermal  charged  particles  are  usually  very  high  along  magnetic  field  lines 
above  the  "sensible"  atmosphere.  Accordingly,  the  thermal-plasma  constituents 
become  redistributed  rapidly  along  the  field  lines  and  tend  to  maintain  the 
potential  the  same  all  along  the  field  line.  This  approximation,  as  discussed 
in  Section  III,  is  good  for  low  and  intermediate-yield  bursts.  For  high-yield 
bursts,  however,  this  assumption  may  overstate  the  debris  at  low  altitudes 
that  drifts  to  higher  L values. 

The  equipotential  field-line  approximation,  however,  greatly  simplifies  the 
computation  of  the  debris  motion.  If  the  field  lines  are  equipotentials, 
the  debris  drifts  outward  to  higher  L values  as  though  the  magnetic  tube  on 
which  the  debris  is  injected  convects  upward  preserving  its  magnetic  flux 
and  dipole-like  configuration.  This  behavior  of  the  plasma  can  be  shown  as 
follows.  Consider  the  debris  to  be  located  initially  in  a magnetic  tube  which 
has  a center  line  that  crosses  the  equator  at  the  geocentric  distance  rQ  (see 
Figure  2).  The  argument  above  requires  that  if  the  field  lines  are 
equipotentials  and  the  electric  field  in  the  tube  is  directed  toward  the  east, 
every  element  of  the  plasma  along  the  tube  will  drift  in  the  local  E X B / , 

direction  at  the  speed  E/B  such  that  it  will  reach  a magnetic  tube  which  has  a 
center  line  that  crosses  the  equator  at  rQ+  Ar0  in  the  time  At.  It  is  therefore 
necessary  to  show  that. 
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= AC 
E/B 


Ar 


E / B 
o o 


(1) 
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where  is  the  displacement  of  the  plasma  element  that  was  located  along  the 
initial  tube  at  the  coordinates  4,  . The  subscript  o refers  to  quantities  at 
the  equator  (0  * 90°).  Conservation  of  magnetic  flux  in  a dipole  field 
requires  that. 


B AC  r sin  0 A$  ■ B At  r 
o o o 


(2) 


or 


AC 


B r 

o o . 

B r sin  0 o 


(3) 


In  Eq.  (2)  A<f>  is  the  azimuthal  width  of  the  tube.  Since  the  field  lines 
are  equipotentials , 

V = e r sin  0 A$  ■ EQrQA$  (4) 


or 


E r 
o o 

r sin  0 


(5) 


Hence,  by  substituting  (3)  and  (S)  in  (1),  Eq.  (1)  is  verified.  Note  that  the 
above  equations  hold  for  any  axisymmetric  magnetic  field. 


15 


b.  Equivalent  Electrical  Circuit 

The  model  of  the  debris  motion  is  further  simplified  by  using  an  equivalent 
electrical  circuit  to  compute  the  excess  charge,  hence  the  electric  field 
across  the  tube.  The  circuit  is  used  principally  to  account  for  the  ionospheric 
current  which  tends  to  neutralize  the  charge  across  the  tube.  It  also  accounts 
for  the  polarization  of  the  thermal  plasma  at  high  altitudes  that  also  reduces 
the  excess  charge.  The  circuit  is  shown  in  Figure  3. 


Figure  3.  Equivalent  electrical  circuit. 


Here,  C is  the  capacitance  of  the  tube  (the  excess  charge  on  the  tube 
boundaries  shown  in  Figure  1 divided  by  the  voltage  across  the  tube);  is 
the  inductance  of  the  current-flow  path  shown  in  Figure  1,  from  the  equator 
to  the  ionosphere;  and  R is  the  effective  resistance  of  the  current  path. 
The  current  Ip  is  the  azimuthal  drift  current  due  to  the  injected  electrons 
and  energetic  ions;  Ic  is  the  current  which  establishes  the  excess  charge 
across  the  tube  (the  current  through  the  capacitor);  and  I^is  the  field- 
aligned  current  which  tends  to  neutralize  the  excess  charge.  The  effect  of 
the  polarization  current  is  taken  into  account  as  discussed  below  by  using 
an  appropriate  dielectric  constant. 
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(1)  Capacitance 


The  capacitance  of  a small  section  of  the  tube  of  length  AS  along  the 
magnetic  field,  "height"  A£  normal  to  B in  the  magnetic  meridian,  and  "width" 
r sin  0 Aip,  is  given  by  the  equation. 


AC  * kk 


AS  AS 

o r sin  9 txp 


where  k is  the  dielectric  constant  of  the  plasma  and  k is  the  permittivity  of 

-12  ° 

free  space  (8.85  x 10  farads/m).  For  the  time  scale  of  the  tube  motion,  the 
effective  dielectric  constant  of  the  plasma  is, 

k - 1 + p/k0B2  (7) 

3 

where  p is  the  density  of  the  plasma  in  kg/m  and  B is  the  magnetic  field 
2 2 

intensity  in  W/m  . Since  P/kQB  » 1, 

k k rs  -8—  . (8) 

° B2 

From  Eq.  (3),  the  ratio  of  the  "height"  to  the  width  of  the  tube 
section  Is 


r sin  6 ’ B sin6 9 ro*P 

when  the  dipole-field  expression  for  r is  used: 

r * r sin29  . 
o 


But,  again  for  a dipole  field, 


B sin69  * B (1+3  cos29)^. 
o 


Hence, 


AS  m 3i  2 — r °— 

r sin  9 A9  (1  + 3 cos  9)*  rQA{p 


17 


■l^J  !,!■■<* 


Since  the  magnetic  bubble  is  nearly  spherical  in  shape,  the  ratio  of  the 
"height"  to  the  width  of  the  section  containing  the  burst  point  (0=9*)  is 
taken  to  be  1 for  the  initial  position  of  the  tube.  Thus, 


Ar 


2 ^ 
(1+3  cos  9*)  . 


(13) 


When  the  tube  moves  to  L values  away  from  its  initial  location,  Lq,  conserva- 
tion of  magnetic  flux  requires  that  ArQ  increase  as  whereas  rQ  Acp  increases 
as  L.  Therefore,  in  general, 


r sin  6 


/ 1+3  cos20*  \ L_ 
'1+3  cos^0  > Lo 


(14) 


By  substituting  (8)  and  (14)  into  (6)  and  summing,  the  initial  capaci- 
tance of  the  tube  is  given  by  the  equatio^ 


C(t=0) 


■E  ^4  1 1 " 3 C-°S'-k- ) 

3 cos  8,  ' 


fk  / il 
B?  W + 


AS 


k* 


(15) 


At  later  times  it  is  necessary  to  take  into  account  the  fact  that  the 
initial  ambient  plasma  in  each  section  of  the  tube  remains  approximately  in 
that  section  while  the  tube  moves.  Hence,  the  ambient  plasma  mass, 

AS^  « AS^/B^  is  a constant  of  the  motion.  (A^,  the  cross  sectional 

area  of  the  tube  section,  is  proportional  to  l/B^).  Hence,  the  general 
expression  for  the  tube  capacitance  is, 

,<A) 

c(t)  -E  ! — - 1 • ~ + 


■f  * 


B,. 


2.  * 

1+3  cos  9^ 


(- 

V1  + 


2 

3 cos  8, 


)'  t 


(16) 


where  p5^  is  the  density  of  the  energetic  ions  in  the  kth  section.  All  of 
K ^ 

the  quantities  in  this  expression  vary  with  time  with  the  exceptions  of  , 

L and  the  invariants  within  the  brackets  designated  by  the  subscript  t=0. 
o 


A model  of  the  ambient  plasma  density  in  the  magnetosphere  with  its 
dependence  on  the  local  time,  season,  and  phase  of  solar  activity  has  been 
developed  to  compute  the  tube  capacitance.  The  model  incorporates  a computer 
program  written  by  Dr.  Y.T.  Chiu  of  the  Aerospace  Corporation  (Ref.  10)  that 
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computes  the  electron  number  density  In  the  altitude  range  90-500  km. 

This  code  includes  the  diural,  annual,  and  solar-activity  cycle  variations. 

The  electron  number  density  is  then  set  equal  to  the  sum  of  the  number  den-  »«; 

sities  of  the  principal  ions  in  that  altitude  range,  i.e.  0^,  N0+,  and  0+, 
and  the  plasma  density  is  computed  by  using  published  data  (e.g.  Johnson, 

Ref.  11)  on  the  fractional  abundances  of  the  ions. 

At  altitudes  above  500  km  extensive  data  on  ion-densities  are  available 
only  from  satellite  measurements  at  altitudes  near  1000  km  and  at  regions 
near  the  equator.  Accordingly,  the  following  procedure  is  used  to  compute 
the  ion  density. 

The  L value  of  the  plasmapause  is  determined  from  the  whistler  data 
(Ref. 12)  on  the  mean  location  of  the  plasmapause  as  a function  of  local  time. 

For  L values  below  the  plasmapause,  the  ion  densities  (0+,  He+,  and  H+)  are 
extrapolated  upward  along  the  magnetic  field  assuming  diffusive  equilibrium 
and  computing  scale  heights  based  on  the  ion  temperatures  at  500  km.  The 
same  method  is  used  for  the  altitude  region  500  km  s h s h(  beyond  the  plasma- 
pause, where  h is  the  "transition"  altitude,  the  altitude  in  the  mid-latitude 

c + 
trough  where  the  principal  ionic  constituent  changes  from  H at  the  lower  lati- 

j. 

tudes  to  0 at  the  higher  latitudes.  The  transition  altitude  is  computed 
from  a model  based  on  the  work  of  Titheridge  (Ref. 13).  At  latitudes  beyond 
the  0+/H+  transition  region,  the  0+  ion  density  is  computed  from  the  scale 
height  as  discussed  above,  but  the  He'*'  and  H*"  densities  are  interpolated 
along  the  magnetic  field  using  the  satellite  measurements  at  the  equator 
and  near  1000  km.  In  the  interpolation  procedure  a Bn  variation  of  the 
plasma  density  along  the  magnetic  field  is  assumed. 

Typical  results  given  by  the  model  are  shown  in  Figures  4 through  7 
for  a variety  of  geographic  positions,  local  times,  seasons,  and  solar 
activity.  Figure  4 shows  the  computed  concentrations  of  O2  , NO  , 0 , He  , 
and  H+  ions,  and  the  total  ion  concentration,  Nj+,  at  altitudes  below  500  km. 

The  results  are  for  the  local  time  1200,  the  month  of  December,  and  the 
geographic  coordinates  of  0°  latitude  and  280°E  longitude.  Figure  5 
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illustrates  the  radial  variation  of  the  total  ion  concentration  at  the  equator 

for  several  combinations  of  the  available  parameters  (local  time,  month, 

geographic  position,  and  solar  activity).  This  figure  clearly  indicates 

+ + + 

the  plasmapause  positions.  Figures  6 and  7 illustrate  the  0 , He  , and  H 
ion  concentrations  in  a geographic  meridian  as  a function  of  magnetic  latitude 
at  altitudes  of  500  and  1000  km,  respectively.  The  high-latitude,  light-ion 
trough  is  visible  in  Figure  7.  Here,  the  light  ions,  H+  and  He+,  are  depleted 
through  their  outward  flow  and  convection  out  of  the  magnetosphere,  and  the  0 
ions  become  dominant  at  1000  km. 

The  tube  capacitance  due  to  the  ambient  plasma  alone  is  shown  as  a 
function  of  L in  Figure  8.  The  results  are  shown  for  local  times  of  noon  and 
midnight,  the  months  of  December  and  June  during  moderate  solar  activity, 
and  for  the  geographic  longitude  of  100°E.  As  expected,  the  capacitance 
rapidly  increases  as  the  lengths  of  the  field  lines  that  thread  plasma  increase 

above  the  100  km  level,  reaches  a local  maximum  for  field  lines  that  lie  in 
much  of  the  F-layer,  then  decreases  above  the  F-layer  as  shorter  sections  of 
the  field  lines  lie  in  the  dense  ionosphere.  At  higher  L values  the  B 
factor  in  the  denominator  of  Equation  (16)  causes  the  capacitance  to  increase 
again  beyond  the  local  minimum  at  1.5  < L < 1.7.  The  sharp  reduction  near 
L=4  is,  of  course,  due  to  plasma  reduction  beyond  the  plasmapause.  In  addition 
to  the  differences  in  the  capacitances  shown  in  the  figure  due  to  the  local 
time  and  season,  the  capacitances  are  also  sensitive,  especially  at  L ^ 2.5, 
to  geographic  longitude.  Near  the  Atlantic  geomagnetic  anomaly,  the  capaci- 
tances for  L < 2 are  lower  than  those  shown  in  Figure  8 by  about  an  order  of 
magnitude . 


(2)  Inductance 

If  the  surfaces  of  the  tube  containing  the  excess  charge  are  regarded 

to  be  rigid  conductors,  the  inductance  of  the  circuit  carrying  the  neutralizing 

current  (See  Figure  1)  is,  approximately. 
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: L-value  for  the  local  times  noon  and  midnight  and 
ig  a geographic  longitude  of  100°E. 


where  nQ  is  the  permeability  of  free  space  (4n  x 10  H/m),  and  SQ  is 

the  length  in  meters  of  the  field  line  from  the  E-region  of  the  ionosphere 
to  the  equator.  For  a dipole  magnetic  field. 


c = _£ 

o 2 


~ ( x Vl  + 3 x2  + jj-  la  C/T  x + Jl  + 3x2)  ] 


where  rQ  is  the  geocentric  distance  to  the  field  line  crossing  of  the  equa- 


tor, and  x * sin  A. 


(^) 


where  h is  the  altitude  of  the  E-region  of  the  ionosphere  («  110  km),rE  is 
the  radius  of  the  earth,  and  L is  the  L value  of  the  field  line. 


(3)  Resistance 

(a)  Undisturbed  Ionosphere 

When  the  tube  has  drifted  to  regions  away  from  the  burst  where  the 
ionosphere  is  not  appreciably  perturbed  by  the  nuclear  explosion,  the 
following  equations  are  used  to  compute  the  resistance  R: 


iv  - , 2 K J 

El+E2 

where  and  E2  are  the  height  integrated  Pederson  and  Hall  conductivities, 
respectively. 
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-20 

where  E^  end  E ^ are  1°  abmohs,  and  o ^ are  in  abmohs/cm,  e = 1.602  x 10  emu, 

n ■ E n, , T is  the  electron  temperature  in  °K,  m.  * ion  mass  in  gm,  m = 
e I i’  e2g  r i e 

9.1058  x 10  gm,  1/  is  the  electron-neutral  collision  frequency, 

is  the  electron-ion  collision  frequency,  is  the  ion-neutral  collision 

frequency,  M is  the  mean  molecular  weight  of  the  ions,  pn  is  the  neutral 

density  in  gm/cm3,  B is  the  magnetic-field  intensity  in  gauss,  u)e  is  the 

electron-cyclotron  frequency,  iu,  is  the  ion-cyclotron  frequency,  W.  is  the 

1 5*5 

molecular  weight  of  the  i th  ion  species,  h,  * 80  x 10  cm,  and  h,  * 400  x 10 

1 _9  x 

cm.  Using  these  values,  Eq.  (20)  must  be  multiplied  by  10  to  obtain  R 
in  ohms. 

The  variables  required  for  the  calculation  of  R are  the  geographic 
position  of  the  foot  of  the  flux  tube  containing  the  burst,  the  local  time, 
the  10  cm  solar  flux,  and  the  temperature  Tg  (this  is  assumed  to  remain 
constant) . 
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Eq.  (20)  gives  values  near  the  equator  that  vary  from  about  0.7  ohms  to  0.06 
ohms  through  the  local  time  sector  0-6  hrs. 
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(b)  Disturbed  Ionosphere 

A nuclear  burst  at  altitudes  within  several  thousand  kilometers  decreases 
this  resistance  by  many  orders  of  magnitude.  Recombination  proceeds  rapidly 
at  low  altitudes,  but  at  the  altitudes  of  the  effective  E-region  of  the 
ionosphere  the  conductivity  remains  extremely  high  for  tens  to  hundreds  of 
seconds.  For  the  equivalent  electrical  circuit  under  consideration  — a series 
resonant  circuit  — the  resistance  does  not  appreciably  influence  the  voltage 
across  the  capacitor,  hence  the  tube  motion,  if. 


R < < 


(33) 


Note  that  the  period  of  the  resonant  circuit  is: 


4 „ e“o  _p_  -,1/2 


t - 2n  C/2  2tt  t -q  m.q  Sq  f — jj  ds  ] 

o B 


(34) 


2 

And  for  p/B  = constant  along  the  field,  Eq.  (34)  becomes, 


t = 4 S / V. 

O 9 A 


where  is  the  Alfven  velocity, 


\ - B/VJJ 


Hence,  the  inequality  (33)  becomes, 


E<<ff  TA  - 1-6  * 1°‘6  ta 


(35) 


(36) 


(37) 


V.  is  about  106  m/sec  in  the  natural  magnetosphere,  but  the  injected  ions  will 
reduce  the  Alfven  velocity  in  the  debris  tube,  perhaps  by  as  much  as  an  order  of 
magnitude.  Therefore,  the  tube  motion  will  be  insensitive  to  the  value  of  R 
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if  R < < ~ .2  0,  i.e. , if  R is  about  0.02  0 or  less. 

This  circumstance  is  fortunate  because  the  resistance  of  the  ionospheric 
path  after  a nuclear  explosion  is  difficult  to  compute  accurately  and  it  changes 
rapidly  in  time.  Eq.  (37)  indicates  that  the  resistance  may  be  set  equal  to 
zero  until  the  foot  of  the  debris  tube  in  the  ionosphere  moves  away  from  the  sub- 
burst point  by  about  1500  km,  where  the  ionosphere  would  not  be  disturbed 
appreciably  and  the  value  of  R given  by  Eq.  (20)  would  apply. 

(c)  Anomalous  Resistivity 

As  mentioned  in  the  introduction,  anomalous  resistivity  may  develop  along 

the  current  path  above  the  ionosphere  if  the  burst  yield  is  sufficiently  high. 

■5 

Kindel  and  Kennel  (Ref.  5)  have  shown  that  field-aligned  currents  of  2 X 10 

c p 

to  2 X 10”D  Amps/m  can  destabilize  ion  cyclotron  waves  within  a wide  range  of 

the  electron  to  ion  temperature  ratio  and  for  critical  drift  velocities  of  the 

current  carriers  that  are  a snail  fraction  ( > O.l)  of  the  electron  thermal 

velocity.  The  resulting  ion  cyclotron  turbulence  would  produce  anomalous 

resistivity  and  cause  high  electrostatic  potential  differences  to  appear  along 

the  magnetic  field  across  the  turbulent  region.  These  predictions  have  been 

verified  recently  in  the  auroral  region  with  data  obtained  with  the  Mr  Force 

S3-3  satellite.  On  repeated  passes  at  high  latitudes,  simultaneous  observations 

"6  2 

are  made  of  field-aligned  currents  of  ~ 10-  Amps/m  , electrostatic  ion  cyclotron 
waves,  plasma  turbulence,  and  electrostatic  potential  differences  of  several  kV 
along  the  magnetic  field  (Refs.  6-8).  in  analyzing  these  data,  Hudson  et  al. 
(Ref.  9)  infer  an  anomalous  resistivity  1\  = 100  ohm-m  compared  with  the  classical 
resistivity  at  the  satellite  altitude  of  1.5  X 10  ohm-ra. 

For  the  MT-yield  test  devices,  the  field  aligned  current  was  of  the  order 
of  10^  Amps/(150  km  X 15  km)  « 5 X 10 Amps/m2.  Hence,  the  essential  criterion 
for  the  instability  is  met.  If  we  assume  the  anomalous  resistivity  to  be  100 
ohm-m  and  to  extend  for  a distance  d « 1000  km  along  the  field  as  observed  in 
the  auroral  zone,  then  the  total  resistance  along  the  field  would  be  R = T|d/A 
a,  100  X 10^m/2  X 109  m2  = .05  ohms.  This  value  of  R is  so  small  that  it  will 
not  appreciably  affect  the  tube  motion,  and  it  is  not  taken  into  account  in  the 
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injection  model.  However,  the  potential  drop  along  the  field  established  by  the 
resistivity  reduces  the  E X B drift  motion  of  the  debris  at  low  altitudes 
relative  to  that  computed  by  the  model.  Uiis  effect  is  discussed  in  Section  III. 


c.  Motion  of  Debris  Tube 

As  discussed  in  Section  II. 2. a.,  the  assumption  that  the  magnetic  field  lines 
are  equipotentials  allows  one  to  determine  the  cross-L  distribution  of  the 
debris  by  computing  the  motion  of  the  magnetic  flux  tube  containing  the  debris. 
Using  the  symbols  defined  in  Sectiai  II-2.a.t  the  tube  moves  in  L at  the  rate. 


Since, 
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V l3 


(38) 

(39) 

(40) 
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where  Bj,  is  the  field  intensity  at  the  equator  (3.12  X 10  J W/m  ),  Eq.  (38) 
may  be  written  as, 
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dt 


rE  BE  ° A 'P 


(42) 


q is  the  net  charge  across  the  tube.  It  is  given  by  the  solution  of  the  circuit 
equation, 


dt 


dt 
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C f dt 
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As  discussed  above,  L , R,  C,  as  well  as  1^,  are  , in  general,  functions  of 


time. 


In  an  attempt  to  verify  a code  (TUBE)  that  was  written  to  compute  the 


tube  motion,  the  values  of  L^,  R,  and  C were  taken  to  be  constants,  with 
R/2Ii£  < (Lf^)  and  the  right  hand  side  of  (43)  was  also  set  equal  to  a 
constant,  F.  The  exact  solution  of  (43)  is  then. 
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Q = [ 1 - ~ — (0!  sin  art  + in  cos  u)t)] 


(44) 


P P 

where  a = R/2L^  and  ou  = 2 /L^C-a  and  the  exact  solution  of  (42)  is 
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The  results  of  this  test  are  given  in  Table  I which  list  the  exact  value  of 
L computed  from  Eq.  (45  ) with  the  numerical  solution  obtained  by  TUBE  with  a 
time  step  of  0.3  sec.  The  parameters  were  R = 1 ft,  Lf  * 4.64  h,  and  C = 10  f. 

Table  1 

Comparison  of  Numerical  and  Exact  Solutions  of  Differential  Equation 


t 

sec. 

L 

exact 

L 

computed 

AL 

error 

0.0 

1.8000 

1.8000 

.0000 

6.0 

1.8096 

1.8096 

.0000 

12.0 

1.8594 

1.8594 

.0000 

24.0 

1.9714 

1.9716 

.0002 

36.0 

2.0058 

2.0060 

.0002 

48.0 

2.1571 

2.1574 

.0003 

60.0 

2.2018 

2.2023 

.0005 
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d.  Azimuthal  Drift  Current 


(1)  Drift  Current  Due  to  Ionized  Debris 
The  azimuthal  drift  current  of  the  hot  ions  (debris  and  hot  air)  at  the 
western  face  of  the  debris  tube  is, 

4 - «>  iN-  *d  (46) 

where  q is  the  mean  charge  of  the  hot  ions,  dN/d^  is  the  azimuthal  gradient 

0 

of  the  energetic  ions,  and  is  the  mean  azimuthal  drift  velocity  (radians/ 
sec)  of  the  hot  ions.  1^  is  in  Amperes  if  q is  in  Coulombs.  The  hot  ions  will 
be  assumed  to  be  uniform  in  every  cross  sectional  area  A,  of  the  tube  normal 
to  the  field.  Hence,  Eq.  (46)  can  be  written  as, 

XD  = 4‘T  / nH  (‘V’  S ^ % ^V’  ds  (*7) 

where  A $ is  the  width  of  the  debris  tube,  n (v,s)  is  the  number  density  of 

H 

the  hot  ions  at  point  s along  the  magnetic  field  (ion  velocity  is  v at  that 
point),  A(s)  is  the  cross  sectional  area  of  the  tube  at  s,  and  4>D  (v,  s)  is 
the  azimuthal-drift  velocity  of  the  hot  ions  at  s.  The  integration  is  over 
the  entire  length  of  the  tube  above  an  altitude  of  200  km. 

If  a net  force  F acts  on  a charged  particle  over  its  gyro  motion,  the  particle 
drifts  with  a velocity  F x B/q  B . 

The  hot -ion  velocity,  as  discussed  in  Section  II. 3,  is  principally  along  the 

2 

field,  hence  the  force  F is  essentially  equal  to  the  centrifugal  force  Mv  / R. 


Thus, 


• - (*£-)  I 

D \ R / qBrcos 


where  M is  the  mean  mass  of  the  hot  ions,  and  R is  the  radius  of  curvature  of 
the  field  line  at  s.  Note  that  the  particle  charge  cancels  when  *D  is  put 
into  Eq.  (47).  Hence,  information  on  the  degree  of  ionization  of  the  hot 
ions  is  not  required. 


For  a dipole  field, 
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Hence, 
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|‘  a v A(s)  f(u)  ds 
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where. 


f(u)  = 


1 - u 

(1  + 3u y 


(52) 


The  mean  mass,  number  density,  and  velocity  of  the  hot  ions  at  s are  computed 
as  described  in  Section  II. 3 from  the  ion  distributions  provided  by  NRL  at 
certain  reference  planes,  normal  to  the  magnetic  field,  above  and  below  the 
burst  point. 

Since  the  ions  in  the  burst-point  region,  between  the  reference  planes,  are 
temporarily  trapped,  their  velocities  are  mainly  perpendicular  to  the  magnetic 
field.  The  net  force  F is  then  1/2  MV  /R,  derived  from  the  product  of  the  ion 
magnetic  moment  and  the  perpendicular  gradient  of  the  magnetic  field.  Hence, 
is  just  one-half  that  given  by  Eq.(48).  Moreover,  since  the  interval 
between  the  reference  planes  is  small,  the  integral  (47)  for  the  hot  ions  in 
that  region,  l.e., 

l£  = — r (1/2  Kv2)  nHaf  (u)  ds  (53) 

rE  h L V 

is  taken  to  be, 
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where  W(t)  is  the  kinetic  energy  of  the  hot  ions  remaining  between  the  reference 
planes  at  time  t,  and  f(u*)  is  given  by  Eqs.  (52)  and  (50)  evaluated  at  the  burst 
point  latitude,  X*. 

(2)  Drift  Current  Due  to  Injected  Electrons 

The  azimuthal  drift  current  of  the  electrons  is  somewhat  more  difficult 
to  determine  analytically.  For  this  case  the  general  expression,  (47),  for 
the  electron  drift  current  at  the  eastern  face  of  the  debris  tube,  where 
cp  = cd2  as  shown  in  Figure  9, 


*2 


Figure  9. 


Schematic  diagram  of  tube  cross  section  illustrating 
computation  of  electron  drift  current  at  4>  = 


may  be  written  as. 


03  t(w) 

(t)  = f dw  T ®p(t_t )g(w  j1;-t *)$_(w)dt ' 

Oo“  D 


( 55) 
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where  q is  now  the  electron  charge. 


Sg(t)  = 6 I’  nf  (s,  t)  A (s)  ds 


nf  is  the  number  density  of  the  fission  fragments, 


. (1  . bW/b//2 


is  the  trapping  fraction  of  the  electrons  that  are  emitted  (isotropically) 
at  s,  and  is  the  field  intensity  at  the  local  atmospheric  cut-off  altitude. 

TOie  factor  6 in  Eq.  ( 56  ) accounts  for  the  number  of  electrons  per  fission 
fragment  that  are  emitted  during  the  decay  of  gross  fission  fragments. 


R (t) 


(1  + t)J 


is  the  beta  decay  rate  of  the  fission  fragments  in  electrons/sec,  normalized 
to  one  electron,  i.e.. 


‘ R(t)  dr,  = 1 

• o 


G (w,  t)  = — 


a-v/v0(t) 

^ctr- 


is  taken  to  be  the  energy  spectrum  of  tine  electrons  in  electrons/MeV.  Here, 
the  characteristic  energy,  wQ(t),  is  a monotoni rally  decreasing  function  of 
time.  The  spectrum  is  normalized  to  one  electron,  i.e., 

30 

f G(w,  t)  dw  = l (6l) 


The  quantity  SeRG/A<|>,  therefore,  gives  the  number  of  trapped  electrons 
that  are  injected  in  the  debris  tube  per  second  per  MeV  per  unit  azimuthal 


angle.  Multiplication  by  q and  by  $Q( v),  the  mean  azimuthal  drift  velocity  of 
electrons  of  energy  w,  gives  the  azimuthal  drift  current  due  to  that  increment 
of  the  electron  source.  The  double  integral,  (55),  sums  the  contributions  to 
the  current  at  t and  cp  = cp2  the  electrons  which  are  injected  at  cp1  :£  cp  £ cp^ 
at  the  earlier  times  t - t* , where 

<P2-<P 

= rrr~  ^ 

(w) 

the  upper  limit  of  the  t'  integration  is 


, t ^ -AJ£— 

«D  (w) 

, t > -A2_ 
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I (w)  -l 


«D  (w) 


Since  the  electron  energy  (=2%  of  fission  yield)  is  much  less  than  the  hot-ion  energy 

the  drift  current  due  to  the  electrons  is  much  less  than  that  of  the  hot  ions. 

Hence,  it  need  not  be  computed  very  accurately,  and  the  following  approximations 
seem  justified, 

S (t  - tO  = S (t)  (64) 

€ 6 


w = <w  (t)  >.  twlMeV 
o o t 


and  ft  (w)  = a w L = 2.44  X 10_J  w L (66) 

In  addition,  by  using  Eq.  (62)  to  change  the  variable  of  integration  from 
t ’ to  cp , Eq,  ( 55)  becomes , 

0.2qSa(t)  « . <P2  V*  \ -1.2 


«p2  - $D  (w)  ts: 


where  cp. 


A cp 

«D(V) 


, t>  -Ai£- 

$D(W) 


The  cp  integration  gives, 
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Furthermore,  since 


*2  " , „ . A co 

STT  = * forw5wi  = ii^- 

$D  (w) 


Eq.  (67)  may  he  written  as, 


qS  (t)  ( W1  / 0 

^ (t)  * v^f-?  - | J'0  !1  - 1 »D  (B)  dw  + 

f e 0 T(l+t  - -^-2-  )"*2-(l+t)",2|  iD(w)  dw  \ (70) 
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+ R (t)  e t ) (71) 


where  t(L)  = - 


V vL) 

is  the  time  for  an  electron  of  energy  wq  to  drift  across  the  debris  tube. 
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DISTRIBUTION  OF  IONIZED  DEBRIS 


I : 


: 

* 

■ i 

li 


As  mentioned  in  the  introduction,  the  distribution  of  the  debris  and  the  i 

accompanying  energetic,  ionized  air  in  the  vicinity  of  the  magnetic  bubble  have 

been  provided  by  NRL.  Initially,  however,  this  ion  distribution  included  only 

the  debris,  and  the  codes  that  computed  the  subsequent  distribution  of  the  ions 

were  based  on  the  analyses  described  below  in  subsections  a.  and  b.  for  the 

debris  alone.  The  modifications  of  the  codes  to  include  the  effects  of  the 

hot  air  are  described  in  subsection  c. 

a.  Initial  Motion  of  Debris  Along  Magnetic  Field 

The  Starfish  debris  data  provided  by  NRL  specify  certain  properties  of  the 
debris  in  reference  planes  perpendicular  to  B above  and  below  the  burst  point. 

For  each  time  step  of  0.1  sec,  0 <_  t >_  10  sec,  the  data  are  given  at  various 
radial  distances,  r^,  from  the  central  field  line  of  the  tube;  the  data  are 
cylindrically  symmetrical  about  this  field  line.  The  properties  which  are 
pertinent  to  the  injection  process  are  the  mass  flux,  dM/da,  and  the  mean 
values  of  the  velocities  parallel  and  perpendicular  to  the  field,  v^  and  Vj^. 

Here  dM/da  is  the  debris  mass  in  grams  per  second  per  cm^  of  the  reference 
plane.  The  reference  planes  are  550  km  above  and  below  the  burst. 

The  data  do  not  include  the  fraction  of  the  debris  that  escapes  through  the 
"neck1'  of  the  magnetic  bubble  at  early  times.  In  a private  communication, 

Dr.  R.  Clark  of  NRL  notified  us  that  only  about  1.5  percent  of  the  debris 
escapes  in  this  manner,  and  that  it  becomes  distributed  over  a radial  distance 
of  about  10  km  at  the  550  km  reference  planes. 

The  data  were  received  on  magnetic  tape,  and  were,  therefore,  easily  processed 
for  inclusion  in  the  model.  First,  the  data  were  averaged  over  5 time  steps 
(0.5  sec).  This  averaging  was  necessary  because  the  data  on  the  tape  were 
quite  noisy.  Examples  of  the  resulting  mass  flux  that  crosses  the  upper 
reference  plane  at  various  time  intervals  are  shown  in  Figures  10  - 12.  Here, 
the  debris  mass  flux,  integrated  over  the  0.5  second  intervals  shown  in  the 
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Figure  1C  Integral  of  debris  mass  flux  from  0 to  0.5  sec  as 
in  upper  reference  plane. 
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tegral  of  debris  mass  flux  from  4.5  to  5 sec.  as 
upper  reference  plane. 


tegral  of  debris  mass  flux  from  9.5  to  10  sec.  as  function  of  radial  distance 
upper  reference  plane. 


figures,  is  plotted  against  the  radial  distance  from  the  central  field  line 
in  the  reference  plane.  Note  that  at  early  times  (0-0.5  sec)  the  flux  is  much 

higher  near  the  central  field  line  and  it  extends  to  smaller  radial  distances 
than  at  later  times.  The  velocities  of  the  debris  do  not  vary  more  than  by 
about  a factor  of  five,  either  as  a function  of  r at  a given  time  or  as  a 
function  of  time  at  r = 0.  Furthermore,  the  v^  components  were  greater  than 
the  components  by  more  than  an  order  of  magnitude. 

The  data  were  then  averaged  to  determine  an  effective  width  of  the  debris 
tube  over  which  the  mass  flux  and  velocity  were  constant  for  each  time  step. 
Since  the  injection  model  utilizes  a tube  cross  section  that  is  square  near 
the  burst  point,  the  distribution  of  the  debris  mass  was  first  determined  as 
a function  of  a rectangular  coordinate  x.  This  distribution  was  obtained 
from  the  integral. 


(73) 


Where  M.  (Jx2+y2' , t„ ) is  the  mass  per  unit  area  crossing  the  reference  plane 

A 1 /r  jr* 

at  the  distance  r = Jx  +y  from  the  central  field  line  in  the  interval  0.5  sec 

and  at  the  time  t . . The  quantity  dx  is  thus  the  debris  mass  that 

1 dx 

crosses  one  of  the  reference  planes  at  x in  the  interval  dx  at  the  time  t^ 
in  the  interval  0.-5  sec. 

The  effective  width,  X+,  of  the  area  in  the  reference  plane  above  the  burst 
and  the  width,  X , of  the  area  in  the  plane  below  the  burst  were  then  determined 
from  the  equation. 


+ 

X"  = 


+ + + 

2 E A M“  x (x£w)±  + 0.0075  x 20  km 


(l  + 0.0075)  M* 


(7*0 
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where  the  superscripts  + and  - refer  to  values  in  the  upper  and  lower  reference 
planes  at  the  time  t^  t .25  secj  (x^w)^  is  the  half -width  at  half -height  of 
the  distribution  given  by  Eq.  (73)  for  the  same  time  interval;  and  is  the 

total  debris  mass  that  traverses  the  upper  reference  plane.  The  mass  is  increased 
by  .7 % to  account  for  the  debris  that  escapes  through  the  neck  of  the  bubble. 
These  quantities  are  given  by  the  equations. 


/ dM"(x, 

V 

l dx 
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max 
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20  4 
E A m; 
i=i 
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(75) 


x = o 


(76) 


(77) 


The  effective  area  was  then  obtained  from  the  equation. 


A = [(X+  + X")/2f . 


(78) 


The  mean  value  of  the  debris  flux  over  this  area  at  t^  is, 


<F“>  = 


2 a m: 


A m 


(79) 


A1 


± 2 
<F*>  is  the  number  of  debris  particles  per  cm  per  second.  The  mean  mass  of 

a debris  particle  has  been  taken  to  be  the  mass  of  an  aluminum  atom 


The  mean  value  of  the  debris  velocity  at  time  t^  is  given  by  the  equation. 
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"Width— at— half -height  of  mass  flux  distribution  along  rectangular 
coordinate  as  a function  of  time.  Curves  a and  b refer  to  the  distri 
buttons  in  the  upper  and  lower  reference  planes  respectively. 


Figure  14.  Mean  parallel  velocity  component  of  debris  in  (a)  upper 
reference  plane  and  (b)  lower  reference  plane  as  a 
function  of  time. 
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The  half-width-at-half  height  of  the  mass  flux  distribution,  x^w,  as  a 
function  of  time  is  shown  in  Figure  13.  The  curve  marked  a in  this  figure  gives 
the  half  width  of  the  distribution  in  the  upper  reference  plane,  and  curve  b 
gives  the  half  width  of  the  distribution  in  the  lower  reference  plane.  The 
average  values  of  the  widths,  weighted  by  the  mass  flux  are  146  km  for  a and 
109  for  b.  The  mean  velocity  components  of  the  debris  along  the  magnetic 
field  at  the  upper  and  lower  reference  planes  are  shown  in  Figure  14  as  a function 

of  time. 


The  fission  fragments  are  assumed  to  be  uniformly  mixed  with  the  debris. 

Hence,  the  flux  of  fission  fragments  can  be  obtained  by  multiplying  the  debris 
flux  by  the  fraction  T],  where: 

Nff 

Tl  = (81) 

D 

N„„  is  the  total  number  of  fission  fragments  released  by  the  burst  (fission 
yield  x 10  ),  and  is  the  total  number  of  debris  particles  (M^  + M^,)  x 

1.015/m^. 


b.  Distribution  of  Debris  in  Magnetosphere 

The  distribution  of  the  debris  in  the  magnetosphere  is  determined  by 
combining  the  cross -L  motion  of  the  debris  tube  with  the  motion  of  the  debris 
along  the  magnetic  field.  For  this  purpose,  as  well  as  for  the  computation 
of  the  capacitance  of  the  tube,  it  is  convenient  to  use  the  grid  shown  in 
Figure  15 . 
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Figure  15.  Illustration  of  partitioning  of  debris  tube  with  surfaces  normal 
to  the  magnetic  field. 


A family  of  curves, 

rk  = Pk  sinV2  Xk 

orthogonal  to  the  local  dipole  magnetic  field  lines  are  used  to  divide  the 

debris  tube  into  sections  of  length  As^.  The  initial  sections  are  defined  by 

a set  of  N + 1 latitudes  { X.  } such  that  X.  > X,  , with  X,  = cos-1 
1/2  ko  ko  K+l,  o k,o 

[(rE+h0)/I'or-EJ  X and  ^W+l  0 = Here>  h0  *s  the  limiting  altitude  due  to 
the  atmosphere.  As  the  tube  moves  outward,  the  boundaries  of  the  sections 
remain  on  the  orthogonal  curves,  (82),  which  are  defined  by  the  N+l  set  of 
polar  intercepts  { pfc} . Since  the  curves  pass  through  the  predetermined  set 
^-Xko^’  the  pol-ar  intercepts  are  given  by  the  equations, 


T 


L 


Lo  ^ \o>  k*N 


*k  = 


(83) 


k = N+l 


A subroutine  DPLATR  has  been  written  to  compute  X^  given  and  L. 


First,  an  approximate  value,  x£,  is  computed  using  the  equations: 
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k | 

tLog10  u^y'2-?  \ + 1-5» 

* < 75 

n/2  - [(l-.25//\)/ 

, \ 

,i+l 


£ 

where  * (L/p^  • 

The  approximate  value,  X^,  is  then  refined  by  using  Newton's  method  with: 

- f^/(df/d X^)^,  1 = 0,1,,., 

f(X*)  = \ COs4  Xk  " sin'Xk 

4 g i i 

-cos  x£  ( 4u^  cos  X£  sin  X^  + l). 


(85) 


where 


and 


(^1 


(86) 


The  iteration  sequence  is  stopped  when  | (AX*+1)/x£  | < e = 10 


v-6 


Knowing  X^  and  L,  other  useful  quantities  along  the  tube  can  be  computed, 
such  as  the  magnetic  field  intensity  at  the  center  of  a section  which  is 
needed  to  compute  the  tube  capacitance  and  the  electron  trapping  fraction. 
This  field  intensity,  e.g. , is  given  by  the  equation, 


48 


J 


tube  in  regions  1 and  3 is  computed  from  the  initial  conditions  of  the'debris 
in  the  reference  planes.  This  computation  is  somewhat  simplified  because,  as 


mentioned  in  the  previous  section,  the  velocities  v^  are  negligible  In  comparison 
with  v,,.  Moreover,  the  specified  values  of  vM  are  so  high  at  the  reference 
planes  that  v ( ( can  be  regarded  to  be  constant  all  along  the  tube.  The  near 
uniformity  of  v ( , can  be  verified  with  the  equation  that  gives  the  variation 
of  v(1  (s)  along  the  tube;  viz., 


| \ 

1 f 


V,  ,2(s)  =•  v,  ,2(s'  ) + V^2(s' ) 


Here,  the  reference  planes  are  located  at  s',  B(s)  and  B(s' ) are  the  field 
intensities  at  s and  s' , g is  the  gravitational  acceleration  at  the  surface 
of  the  earth,  and  r and  r*  are  the  geocentric  distances  to  s and  s’.  This 
equation  reveals  that  a debris  particle  with  the  lowest  specified  velocity 
components  will  loose  only  a few  percent  of  its  initial  v( , component  as  it 
travels  along  the  limiting  field  line  of  the  trapping  region,  from  near  the 
surface  of  the  earth  to  the  equator.  Accordingly,  by  using  the  equations 
for  the  continuity  of  the  debris  particles  and  the  conservation  of  magnetic 
flux  (BA  = const.)  along  the  tube,  the  mean  value  of  the  debris  number  density 
in  the  tube  sections  of  regions  1 and  3 is  given  by  the  equations: 


F(s',t') 

v(t') 


+ F(s*,t*  + At*) 
v(t'  + fit*) 


(89) 
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where  F is  the  debris  flux  in  particles/om  , and 


(90) 


The  latter  equation  follows  from  the  definitions. 


s(t)  - s'(t')  = v(t')  (t-t') 


(91) 


and 

s(t)  + As(t)  - s'(t')  = v(t'  + At')  (t  - t'  - At')  (92) 

In  region  2,  between  the  reference  planes,  the  number  density  is, 

yt) 

^(t)  > = 

t.  , 

where  NR(t)  = Nf  - E J F*  (t')  A±  (t')  dt'  (94) 

+,-  o 

and  V(t)  is  the  volume  of  the  region  at  the  time  t.  is  the  total  number 

of  fission  fragments  released  by  the  burst,  and  the  integrals  in  (9*0  give 
the  number  of  fission  fragments  which  have  left  the  region  at  the  time  t 
through  the  upper  and  lower  reference  planes. 

A subroutine,  DEBATL,  has  been  written  to  compute,  for  each  time  step  of 
the  tube  motion,  the  distribution  of  the  debris,  as  described  above,  and 
other  quantities  that  are  dependent  on  the  debris  distribution;  namely,  the 
tube  capacitance  due  to  the  debris  loading,  the  integral  (*»-2)  for  the 
computation  of  the  drift  current  due  to  the  debris,  and  the  integral  Sg(t) 
for  the  computation  of  the  drift  current  due  to  the  trapped  electrons. 
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c. 


Code  Modifications  for  Inclusion  of  Hot  Ionized  Air 


r 


The  new  data  supplied  by  NRL  include  the  hot  ionized  air  that  expands  into  the 
debris  tube  at  early  times.  The  data  that  are  now  listed  against  time,  t,  and 
radius,  r,  in  the  reference  planes  are: 

Qjj(t,  r),  the  mass  flux  (g/cnr)  of  hot  ions  (debris  and  ionized  air) 
integrated  over  the  time  step, 

Vn(t,  r),  the  velocity  of  the  hot  ions  parallel  to  the  magnetic  field, 

f(t,  r) , the  ratio  of  the  debris  mass  to  the  hot-ion  mass,  and 

M(t),  the  mean  mass  (A.M.U.)  of  the  ions  averaged  over  the  reference  plane, 
at  the  time  t. 

The  effective  area  of  the  debris  tube  and  the  values  of  ion  velocity,  v(t), 
and  hot-ion  flux,  F(t),  averaged  over  the  reference  planes  at  the  specified 
times  t,  are  computed  as  discussed  in  subsection  a.  above,  except  that 
(^(t,  r)  is  used  instead  of  MA(t,  r)  in  equations  (73),  (76)  and  (80),  and 
the  mean  ion  mass,  M(t)/(Avogadro's  Number),  is  used  instead  of  mAi  in  equation 
(79).  Examples  of  the  results  for  Standard  Spartan  bursts  at  altitudes  of 
200,  300,  400  and  600  km,  all  at  L*3,  are  shown  in  Figures  16  - 18.  The 
input  data  were  the  same  in  the  upper  and  lower  reference  planes;  therefore, 
the  results  apply  for  either  plane.  Figure  16  shows  the  sector  width,  A<J> , of 
the  debris  tube  as  a function  of  the  altitude  of  the  burst.  The  sector  width, 
of  course,  is  determined  from  the  square  root  of  the  effective  area  of  the 
debris  tube  at  the  burst  altitude.  In  Figure  16  is  plotted  on  semi-log 
paper  to  emphasize  the  different  scaling  regime's.  For  bursts  at  altitudes 
less  than  about  350  km,  the  results  indicate  that  the  lateral  expansion  of 
the  debris  is  limited  by  the  atmospheric  density  (the  slope  of  the  straight 
line  is  consistent  with  the  scale  height  of  the  atmosphere  in  the  range  200- 
350  km).  At  altitudes  above  about  400  km  the  debris  expansion  is  evidently 
limited  by  the  pressure  of  the  magnetic  field.  The  velocities  of  the  hot  ions 
along  the  magnetic  field  as  a function  of  time  are  shown  in  Figure  17.  At 


times  less  than  about  0.6  sec,  the  ion  velocities  are  about  the  same  for  all 
of  the  burst  altitudes.  At  such  early  times  the  ions  consist  principally  of 
debris  particles  alone.  At  later  times,  the  velocities  of  the  ions  increase 
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Figure  17 • Velocity  of  ions  along  magnetic  field,  in  ref  erence  plane,  versus 
time  for  Standard  Spartan  bursts  at  altitudes  of  2'0.  300,  400, 
and  600  km  (L  - 3) . 
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as  the  burst  altitude  increases.  The  ion  fluxes  as  a function  of  time  are 
shown  in  Figure  18.  In  this  figure  the  ion  flux  for  the  200-km  burst  is  given 
directly  by  the  ordinate.  The  flux  for  the  other  bursts  is  obtained  by 
dividing  the  ordinate  by  the  factors  given  in  parentheses  at  the  curves.  At 
times  beyond  a few  seconds,  the  ion  flux  consists  mainly  of  air  ions. 


In  the  injection  model,  it  is  necessary  to  compute  the  number  flux  of  fission 
fragments  Fff(t)  averaged  over  the  reference  plane  at  time  t.  This  was  done 
by  computing,  as  a function  of  time,  n(t),  the  ratio  of  the  number  of  fission 
fragments  crossing  the  reference  plane  during  the  time-step  interval  to  the 
corresponding  number  of  hot  ions.  This  is  given  by  the  equation, 


n(t) 


Yf  x 1026 
MT(Debris)/M^j 


‘ < f(t) 


> 


R.P . 


. M(t) 

Mai 


(95) 


where  Yf  is  the  fission  yield  in  MT,  Mj,(Debris)  is  the  total  mass  of  the 
debris,  and  < f(t)  >51. p#  is  the  mass  ratio  f(t,  r)  averaged  over  the  reference 
plane ; i . e . , 


fn  f(t)QM(t,  r)  2wrdr 
< f <t)  >R.P  • - /oo -= 

J0  QM(t’  r)  2lrrdr 


(96) 


Note  that  the  term  in  brackets  in  (95)  is  the  total  number  of  fission  fragments 
divided  by  the  total  number  of  debris  (this  ratio  is  assumed  to  remain  constant 
in  any  sample  of  debris).  Recall  that  Maj,  the  mass  of  aluminum,  is  the  mean 
mass  of  the  debris.  This  mass  cancels  in  (95).  The  factor  outside  the 
bracket  gives  the  ratio  at  time  t of  the  number  of  debris  ions  crossing  the 
reference  plane  during  the  time-step  interval  to  the  corresponding  number  of 
the  hot  ions.  The  sum  of  the  numerator  of  (96)  over  the  time  steps  gives  Rj.(Debris). 
The  fission-fragment  number  flux  is  given  by  the  equation. 


n(t)  L QM(t’  r)  2irrdr 

F»(t>  ’ mftTt 


(97) 


where  A^  is  the  effective  cross  sectional  area  of  the  debris  tube  in  the  reference 
plane,  and  At  is  the  time  step.  An  example  of  the  fission-fragment  flux  is  shown 
in  Figure  19,  where  Fff(t)  and  F(t),  for  the  Argus-3  burst,  are  plotted  against 
time. 
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Figure  18. 


Hot-ions  flu:c  at  reference  planes  vs.  time  for  Standard-Spartan 
bursts  at  altitudes  of  200,  300,  400,  and  600  kn. 




According  to  Dr.  R.  Clark  of  NRL,  the  accuracy  of  the  data  on  the  hot  ions 
diminishes  rapidly  at  times  greater  than  a few  seconds.  This  inaccuracy  in 
the  available  data,  as  well  as  the  incompleteness  of  the  data  (t<  10  sec) , 
introduces  an  uncertainty  in  the  tubo  motion.  The  sensitivity  of  the  tube 
motion  to  different  ion  distributions  has  not  been  fully  assessed.  However, 
from  limited  tests  that  have  been  made  with  and  without  the  hot-air  ions,  it 
appears  that  the  error  in  the  tube  motion  may  not  be  severe. 


The  main  program  DEBCON  and  the  subroutine  DEBATL  have  also  been  modified  to 
include  the  effect  of  the  hot  ionized  air.  DEBCON  now  allows  input  data  that 
specifies  Mj.(t)  (the  mean  mass  M(t)  in  A.M.U.)  and  n(t)  which  vary  with  time. 
Previously,  when  the  NRL  data  included  the  debris  alone,  the  mean  mass  was  a 
constant,  = 27  A.M.U. , and  n was  a constant. 

As  discussed  in  the  previous  section,  the  subroutine  DEBATL  computes  the 
capacitance  of  the  tube,  the  azimuthal  drift  current,  and  the  integral  of  the 
number  density  of  the  fission  fragments  in  the  tube.  The  density  of  the  hot 
plasma  (debris  and  hot  air)  in  the  tube  is  required  to  compute  the  capacitance 
and  the  azimuthal  drift  current.  The  f ission-f ragmei.i  number  density  is 
computed  from  Fff(t)  that  is  specified  in  the  reference  planes.  The  density 
of  the  hot  plasma  at  a point  along  the  tube  is  then  computed  from  the  number 
density  of  the  fission  fragments  at  that  point  by  utilizing  the  values  of 
Mj(t)  and  n(t)  at  the  earlier  time  t when  the  fission  fragments,  moving  the 
velocity  vjj(t),  crossed  the  reference  plane. 
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A.  EARLY-TIME  REDISTRIBUTION  OF  ELECTRONS 


a.  Redistribution  in  Convecting  Debris  Tube 

The  redistribution  of  electrons  in  the  debris  tube  is  appreciable  when  the 
electrostatic  potential  V across  the  tube  is  high;  i.e.,  when  the  tube  is 
moving  outward  rapidly  (L  large).  The  situation  is  depicted  in  the  sketch 
(Figure  20). 

The  electric  field  E in  the  tube  causes  the  electrons  to  drift  outward  at 
the  same  rate  as  the  tube.  However,  the  energetic  electrons  also  drift  toward 
the  east  owing  to  the  configuration  of  the  magnetic  field.  Hence,  in  a 
stationary  frame  of  reference,  electrons  that  are  injected  at  the  western 
edge  of  the  tube,  4>w,  will  follow  a trajectory  such  as  that  shown  by  the 
broken-line  curve  in  the  sketch. 
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Figure  20  Illustration  of  displacement  of  electron  in  L due  to  outward 
drift  motion  of  debris  tube. 


The  trajectory  is  given  by  the  differential  equation, 

dL  = d£  (98) 

L ip 

where  the  dots  indicate  time  derivatives;  ip  * ip  (w,L,u0)  is  the  mean 
azimuthal  drift  velocity  of  an  electron  of  kinetic  energy  w and  equatorial 
pitch-angle  cosine  p0  at  L. 


The  electron  flux  j at  L»q>E, C is  related  to  the  flux  j'  originating  in  the 
tube  at  L'tt'tt',  through  Liouville's  equation  for  relativistic  electrons 
(p2  - w2  + 2m  c2w).  Hence, 


j(L,w,p0,*E,t)  _ j'a'.w’.Pp’^'.t') 
w(w+2mQc2)  w'(w'+2m0c2) 


where  in  c2  is  the  electron  rest  mass  energy,  the  pointt  L,$_  and  L',4>’  are 
0 & 

connected  by  the  trajectory  given  by  Equation  (98), 
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«’  = w + TiJ  *(v 


i (/  is  the  potential  energy,  qV,  in  the  tube  when  the  electron  is  at  $,  i.e., 
at  the  time. 


-I  ? 


(101) 


b yd 


and  viq  and  y ' are  related  through  the  conservation  of  the  first  and  second 
adiabatic  invariants. 


Now,  since  the  fission  fragments  emit  electrons  continuously,  the  flux  at 

L,(f„,t  is  equal  to  the  sum  of  the  fluxes  of  the  electrons  which  are  emitted 

at  (j>  at  time  t plus  those  which  were  emitted  at  earlier  times  along  the 
b 

trajectory  shown  in  Figure  20  and  drifted  to  <(>E,L  at  the  time  t.  Thus,  using 
(99),  the  flux  may  be  written  as, 


j(L,w,y  <)»E,t)  - f E w(w+2m  C2)  dj* 


1 w$(W((,+2moc2)  d4 


d<j>  (102) 


where  the  subscript  <j>  denotes  source  parameters,  and  the  limit 


■■  Maximum  of  | $ at  t = 0 

/ $ where  y = critical  value 


(103) 


But,  as  discussed  in  our  earlier  work  (Reference  1), 


V * R(t<)>)e  f. d8<>  (104) 

dS  ' ' * * * w (tA)  j 2ny  t. 
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where  R(t  ) is  the  emission  rate  at  t normalized  to  one,  i.e.. 


\\ 

1 1 

[-1 

I i 


I' 


R<y 


0.2 


<1+V 


1.2 


(105) 


w^t^)  is  the  characteristic  energy  of  the  source  spectrum  at  t^;  n(s^)  is 
6 times  the  number  density  of  the  fission  fragments  at  the  point  s^  along 
the  debris  tube  (6  is  the  number  of  electrons  released  per  fission  fragment); 
is  related  to  a0<(!  through  the  conservation  of  the  magnetic  moment,  viz.. 


« 


ii 


y = i 


B(s,  <}>) 

Bo<V 


(l-p 


(106) 


where  B(s,<(>)  and  BQ(L(j))  are  the  magnetic  field  intensities  at  s and  at  the 
equator,  respectively  in  the  meridian  <j>,  and  t^  = .L^)  is  the  electron 

bounce  time.  The  integral  is  over  the  length  of  the  debris  tube  between  the 
mirror  points 


! 


it 


Bm(s,<0)  - (107) 

By  substituting  (103)  into  (102)  and  noting  that  the  source  is  uniform  in  <t>. 
Equation  (102)  becomes, 


dj  (L,w,p0,t,<|>E) 
dt 


w(w+2m0c 


2ttA4> 


!)  / p>(s»t|j,)R(t^)Exp  [_-w((|/wo(t()))J  dsd(j> 

c2)  [l-  (1-U02)J  tbw 


(tj 


(108) 


fi  where  A<j>  * $E-<t,W’  an<*  as  indicated  above,  t^  = t (<J>) , 

t (<J>E)  - t,  w^  - w(<j>),  w(<J>E)  * w,  = L(tj>) , L(<J>E)  ■=  L. 

Il 

j Here,  the  mirror  points  of  the  electrons  at  the  different  <p  values  are  determined 

I from  the  following  equation,  which  conserves  the  first  two  adiabatic  invariants 

1 of  the  electron  motion: 
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(109) 


B (S,<j>) 
m 

B (S,* 
m a 


,<p)  pL(4>£)  ' 

Tk)  = h (♦)  . 


where 


l+.853950y-1.35395Qy 

l-.463488y5'4 

„ 1/2 


y = <i-V> 


(111) 


The  flux  given  by  Eq.  (108),  as  well  as  integrations  of  the  flux  over  time, 
equatorial  pitch  angle,  and  energy,  are  computed  by  the  code  DSTRIB  and  as- 
sociated subroutines  described  in  Section  II. 6.  The  integration  of  (108) 
over  the  total  time  required  for  the  debris  tube  to  traverse  the  field  line 
at  L,  i.e.,  over  the  full  time  available  for  the  electrons  to  drift  out  of  the 

tube  onto  L,  times  the  dilution  factor  A4> /2tt  gives  the  equatorial  directional 
-2  -1  -1 

flux  (cm  ’sec  *sr  -MeV)  of  the  electrons  after  they  become  spread  uniformly 
in  longitude.  The  additional  integrations  over  pitch-angle  and  energy  give  the 
equatorial  omnidirectional  flux  with  energies  above  certain  prescribed  minimum 
energies.  Of  course,  knowing  the  directional  flux  at  the  equator,  the  direc- 
tional and  omnidirectional  fluxes  at  any  B value  off  the  equator  can  be  com- 
puted easily. 

The  flux  obtained  from  the  integration  of  (108)  over  time  and  pitch  angle  at 
various  L values  for  the  Starfish  burst  is  shown  in  Figure  21.  Note  that  the 
spectra  of  the  electrons  are  quite  hard  at  low  L values  — more  so  than  the 
equilibrium  fission  beta-decay  spectrum  — and  they  rapidly  become  softer  toward 
higher  L values.  At  large  L values  the  redistribution  integral  includes  most 
or  all  of  the  debris  tube  azimuthal  width.  However,  near  the  beginning  of  the 
debris-tube  motion,  at  the  lower  L values,  the  drift  velocities  of  the  lower- 
energy  electrons  are  so  small  that  the  effective  width  of  the  source  is  reduced 
(<(>^  > (ji^) . This  effect  accounts  for  the  fall»off  of  the  flux  toward  lower 
energies.  The  comparison  of  these  spectra  with  the  measured  spectra  of  the  Star 
fish  electrons  is  discussed  in  Section  III. 

The  flux  at  ><j>E  during  the  early-time  expansion  of  the  trapped  electrons 
around  the  earth  can  be  computed  from  the  flux  given  by  (108)  at  by  using 
Liouville's  theorem,  as  was  done  in  Ref.  1.  However,  the  computation  is  now 
somewhat  more  difficult  because  the  energy  spectrum  of  the  redistributed  elec- 
trons does  not  have  the  simple  exponential  form  that  was  used  previously  to 
perform  the  pertinent  energy  integration  analytically.  A discussion  of  the 
computation  of  this  "drift-diluted"  flux  is  given  below  in  Section  II. 4. 
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Figure  21.  Energy  spectra  of  redistributed  Starfish-injected  electrons  at 
various  L values. 
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b.  Redistribution  East  of  the  Debris  Tube 


t 


/ 


i 

a 


Even  after  the  electrons  drift  out  of  the  tube,  they  will  continue  to  drift 
outward  in  L owing  to  an  eastward  electric  field  which  exists  in  the  region 
of  the  azimuthally-bunched  electrons.  This  electric  field  develops  because 
excess  charge  appears  on  field  lines  due  to  the  azimuthal  drift  of  the  elec- 
trons that  cannot  be  neutralized  immediately.  This  charge  can  only  be  neutral- 
ized by  currents  that  flow  along  the  magnetic  field  lines  containing  the  ex- 
cess charge  and  across  magnetic  field  lines  in  the  ionosphere.  Such  neutrali- 
zation cannot  occur  in  times  less  than  about  the  time,  T^,  required  for  an 
Alfv€n  wave  to  traverse  the  length  of  a field  line.  The  excess  charge  is  roughly 
proportional  to  the  rate  of  change  of  the  number  of  trapped  betas  with  respect 
to  , dN / d<p , times  the  product  of  T and  the  mean  azimuthal  drift  rate  of  the 

A 

electrons  at  <$> . The  charge  is  negative  at  the  leading  edge  of  the  electron 
bunch  and  positive  at  the  trailing  edge.  The  excess  charge  is  positive  at  the 
trailing  edge  because  the  betas  drift  away  from  the  field  lines  on  which  their 
charges  had  been  neutralized. 

The  analysis  of  the  redistribution  of  the  electrons  due  to  this  cause  was  car- 
ried out  only  until  preliminary  results  indicated  that  it  was  not  as  large  as 
that  due  to  the  convection  of  the  debris  tube.  However,  the  process  warrants 
further  effort,  especially  since  the  Telstar  satellite  data  indicate  that  the 
electrons  injected  by  the  high-yield  tests  were  transported  to  higher  L values 
than  can  be  attributed  to  the  convecting  debris  tube  alone. 

5.  DEVELOPMENT  OF  THE  DRIFT  DILUTED  FLUX 

In  the  present  injection  model  the  electrons  are  presumed  to  be  released  when 
they  reach  the  eastern  edge  of  the  debris  tube.  The  source  is  therefore  ap- 
proximately a delta  function  in  (magnetic)  longitude  coincident  with  the  east- 
ern edge.  This  is  in  contrast  with  the  previous  injection  model  where  the  source 
was  spread  out  over  the  entire  tube.  The  principal  drawback  in  that  model  was 
that  the  source  was  assumed  uniform  in  longitude,  with  no  consideration  given 
to  how  the  electrons  actually  leave  the  tube.  Now  the  motion  of  electrons 
within  the  tube  can  be  predicted  by  the  "redistribution"  model.  The  motion  after 
passing  the  eastern  edge  of  the  tube  is  entirely  determined  by  uniform  drift 
motion,  in  conjunction  with  the  loss  processes. 


t There  are  two  principal  methods  for  treating  the  drift  dilution  of  a longitudi- 

nally nonuniform  distribution  of  trapped  electrons.  The  more  straightforward 
method  is  to  follow  bunches  of  electrons  with  a finite  energy  spread.  The  longi- 
tudinal extent  of  a bunch  then  increases  because  of  the  energy  dispersion  of 
the  angular  drift  velocity  4>p,  thus 

d<L 

% (T-t)  Iw-  w (112) 

where  t is  the  injection  time,  T is  the  observation  time,  and  w is  the  kinetic 
energy  of  the  electrons.  A distribution  over  the  longitude  segment,  A4>  can  then 
be  assumed,  and  the  contributions  to  the  total  distribution  added  up.  The 
electrons  may  be  assumed  uniform  over  A<J> ; or,  alternatively,  a Gaussian  dis- 
tribution may  be  used  to  ensure  a smooth  distribution.  This  approach  was 
rejected  for  the  original  Specter  trapping  model  for  two  reasons.  First,  the 


simplest  choices  for  the  distribution  over  A*  are  quite  non-physical,  and  could 
lead  to  serious  errors  if  a bunch  of  electrons  is  followed  for  a complete  cir- 
cuit of  the  earth.  The  second  objection  is  that  following  bunches  of  electrons 
for  1,  2,  3,  4,  etc.  circuits  of  the  earth  can  lead  to  bookkeeping  problems. 

Both  problems  can  be  minimized  by  reassembling  the  complete  energy  and  longi- 

j ; 

!tude  distribution  at  appropriate  time  intervals;  but  their  implications 
remain  very  serious  if  a broad  spectral  range  is  to  be  treated  properly. 

The  second  approach,  which  was  adopted  for  Specter,  is  to  employ  Liouville’s 

, theorem  to  follow  the  evolving  distribution.  Lieuville's  theorem  states  that 

ii 

the  phase  space  density  remains  constant  throughout  the  drift  motion.  Thus  a 
delta  function  distribution  in  energy,  longitude,  and  time  remains  forever  a 
delta  function.  The  longitudinal  spreading  is  accounted  for  by  adding  the  con- 
j tributions  of  an  infinite  number  of  delta  functions.  The  distribution  at  any 

arbitrary  time  is  then  given  by  a single  time  integration  over  the  source  (Ref.  1) . 
' An  apparent  disadvantage  is  that  the  differential  energy  spectrum  at  early  times 

may  have  a very  irregular  form.  This  leads  to  computational  difficulties;  but  the 
■ irregularities  are  actually  real  — if  they  do  not  appear  in  treatments  based  on 

| the  alternate  method,  that  is  a consequence  of  artificial  smoothing  of  the  dis- 

1 tribution.  It  is  worthwhile  to  ask  here  whether  the  Liouville  method  can  be 

altered  or  improved  to  remove  the  computational  difficulties. 
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Consider  a delta  function  (in  time  and  longitude)  source  of  injected  electrons. 
At  any  longitude,  electrons  of  a given  energy  will  pass  any  point  at  regular 
intervals,  corresponding  to  1,  2,  3,  etc.  circuits  of  the  earth.  The  instantan- 
eous energy  spectrum  (see  Figure  22_;  therefore  consists  of  a series  of  delta 
functions,  spaced  at  intervals 

AW  = — (113) 

C d$D 

The  spacing  decreases  with  time  so  that  eventually  the  distribution  appears  to 
be  uniformly  filled,  with  nearly  infinitesimal  gaps.  Obviously,  if  one  wishes 
a smoothened  representation  of  the  spectrum,  one  must  sample  at  intervals  no 
closer  than  the  spacing  given  by  Eq . (113).  Further  smoothing  is  possible  if  it 
is  recognized  that  any  observation  of  a particle  flux  must  take  place  over 
finite  intervals  of  time  and  longitude.  So  a certain  amount  of  time  and  longi- 
tude averaging  — consistent  with  the  duration  of  the  observation  — is  appro- 
priate, and  in  some  cases,  necessary. 

The  energy  spectra  of  the  redistributed  electrons  are  nearly  exponential  and 
can  be  represented  well  by  a sum  of  exponentials: 

j“  c0+ci  exp(-w/Wp)+C2  exp(-2w/wo)+C^  exp(-3w/wQ)+. . . (114) 

(This  is  simply  a Taylor  series  in  exp  (-w/w^).)  Figure  23  demonstrates  how 

well  the  computed  spectra  can  be  fitted  by  as  few  as  3 terms  of  the  form  Cnexp 

(-w/w  ),  with  w, ,wOJ  and  w„  not  commensurate.  A fall-off  at  low  energies  is 
n 1 z J 

accounted  for  by  a subtractive  terms,  — C^exp  (-w/w^).  A straightforward 
prescription  for  constructing  the  spectrum  (114)  is  to  select  a wq  that  gives 
the  correct  assymptotic  behavior  at  large  w.  Then  a least-squares  method  can 
be  used  to  compute  as  many  C^'s  as  desired;  C0is  zero. 

The  contribution  of  the  source  at  time  t to  the  total  flux  of  electrons  above 
an  energy  w^  is  (Ref.  1) 

•00 

Aj  = — EC  E 5 [$-2*m4n(T-t)]exp(-nw/w  ) (115) 

4 w n u o 

o n m 

where  5 is  a Dirac  delta  function  and  m is  the  number  of  completed  circuits 


Figure  22.  Energy  spectra  of  drifting  electrons  at  arbitrary 
instants  of  time,  t^  and  t£  tj_.  The  source  was 
assumed  very  narrow  in  longitude. 
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Figure  23.  Energy  spectra  of  the  injected  electrons  at  several 
L-shells,  as  predicted  by  the  new  injection  code. 

The  dashed  curves  are  approximate  representations  by 
sums  of  exponential  functions. 
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The  drift  velocity  is 


♦ = n . .. 

h 1 wTm" dZ~ 
o 


(116) 


where  is  a function  of  L and  pitch  angle.  To  simplify  the  above  integral. 


as  in  the  treatment  of  Ref.  1,  let 


4>d  % DE 


(117) 


where  the  numerical  value  of  D is  somewhat  different  from  D.  and  E is  the 

2 » 1 
total  energy,  w+m  c , of  the  electron.  This  approximation  is  valid  at  higher 


energies;  at  low  energies,  inclusion  of  the  neglected  factor  leads  to  more 
complicated  results.  In  the  original  model  was  also  assumed  independen 
pitch  angle;  this  simplification  has  physical  implications  (in  contrast  to 
purely  numerical  approximations)  but  it  is  not  really  necessary.  Integra- 
tion of  Eq.  (115)  now  gives  (Ref.  1),  for  a narrow  source  of  width  A$: 

T 


i<Ei> 


* idt  Z K 
= | n n 

a J 


nAj> 


DwQ(T-t) 


x exp{-n($+2irm^)  /Dwq  (T-t)  } / [ l-exp{  -2irn/Dwo  (T-t)  } ] 


(118) 


where  m^  is  an  integer 


DE. (T-t)-$  DE_(T-t)-$ 

1 — i V— ^ +1 


(119) 


2ir 
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The  spectrum  and  source  strength  are  included  in  the  coefficients  K^. 


In  the  existing  Specter  model  the  integration  (118)  is  performed  for  each  point 
E,$,t  in  a grid,  and  is  subsequently  Interpolated  to  obtain  fluxes  at  inter- 
mediate points.  This  procedure  has  several  disadvantages  when  one  has  to  re- 
construct the  differential  energy  spectrum.  The  need  for  interpolations  al- 
ways introduces  numerical  errors.  These  errors  are  compounded  when  the  dif- 
ferential spectrum  at  E is  found  by  taking  the  difference 


% (j(E+AE)-j(E))/AE 


(120) 
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Consider  what  happens  as  a pulse  of  electrons  arrives  at  a point  after  com- 
pleting an  integral  number  of  circuits  of  the  earth.  Figure  24  shows  the  peak 


Tire 


Figure  24.  The  variation  of  integral  fluxes  as  a bunch  of  drifting  electrons 
arrives  at  the  observation  point.  The  figure  is  exaggerated  to  emphasize  how 
the  difference  between  two  integral  fluxes  changes  with  time. 

of  the  pulse  arriving  later  for  the  lower  energy.  Immediately  before  the  pulse, 
however,  there  is  a decrease  in  the  differential  flux  as  the  two  integral  flux 
curves  approach  each  other.  It  is  easy  to  imagine  situations  where  slight 
errors  in  the  interpolation  procedures  could  alter  the  relative  positions  of 
the  integral  flux  curves,  thereby  accentuating  the  irregularities  in  the  dif- 
ferential fluxes. 

One  way  of  dealing  with  spectral  irregularities  is  to  use  relatively  large  energy 
intervals  and  average  over  time,  as  suggested  above.  However,  a more  elegant 
alternative  is  available  that  should  completely  eliminate  the  interpolation 
difficulties.  That  is  to  construct  a formula  that  gives  the  fluxes  as  algebraic 
functions  of  $,  T,  and  E^.  Eq.  (118)  contains  only  one  integration,  and  be- 
comes algebraic  for  a source  that  is  a delta  function  of  time,  t.  The  integra- 
tion has  been  performed  numerically  by  a finite  differencing  method.  The  source 
could  be  replaced  by  a series  of  delta  functions  at  times  t^,  t£»  t^,...t^,  etc., 
with  accuracies  comparable  to  the  present  method.  The  physical  consequences  of 
this  approximation  are  negligible  if  the  spacing  of  the  time  points  is  suf- 
ficiently fine.  The  integral  flux  is  then 


J(V^UKn,k  srrtt-t,)' 

ok  k 


x exp{-n(*+2irmlk)/IX#ok(T-tk)}/[l-exp{-2irn/Dwok(T-tk)  ] 


(121) 


DE  (T-t  )-4>  DE-CT-t.)-* 

-S'  1»1<  - -2  +1 


(122) 


where  WQk  is  the  value  of  wq  appropriate  to  tk;  and  k is  the  contribution 
to  the  n’th  exponential  at  time  tk*  The  differential  flux  at  is 

4i<V  ... 

4 Bj  kn,k  SST  (T-t  > e*P{-»'>+2»">lt)/I)“ot<T-tk)) 


x [exp{+  TrnAm/Dwok(T-tk)}-exp{-7rnAm/Dwok(T-tk)}]/[l-exp{-2Trn/Dwok(T-tk)}]  (123) 


Am  * m1(E1+AE/2)-m1(E1-AE/2) 


(124) 


The  condition  (113)  on  E arises  naturally  in  the  definition  of  m;  if  at 
Ej+  E/w  and  at  E^-  E/w  are  identical,  the  differential  flux  is  zero. 

Eqs.  (121)  and  (123)  give  accurate  results,  are  easy  to  compute,  and  require 
minimal  data  storage.  The  source  might  be  well  represented  by  20  to  50  t^'s 
and  their  associated  w^'s.  Each  tk  might  be  associated  with  a spectrum  that 
can  be  represented  by  10  to  20  exponentails . Thus,  at  each  pitch  angle  or  mir- 
ror field,  the  flux  could  be  adequately  represented  by  a matrix  of  200  to  1000 

K , 's. 
n,k 


6. 


DESCRIPTION  OF  COMPUTER  CODES 


/ 


The  incorporation  of  the  new  tube  motion  debris  model  into  the  SPECTER  codes 
has  resulted  in  the  development  of  3 new  main  programs  AMBCAP,  DEBCON  and  DSTRIB. 

These  main  programs  utilize  some  21  new  subprograms  ACATL,  ARANGD,  BOUNCE, 

BSFRMP,  CONMAG,  DEBATL,  DPLATR,  DTEQU , EXPINT.  IONCON,  LESSEN,  LINTRP,  MAGDEN, 

MOVEL,  PHIDOT,  PLACE,  PPLOC  , SINMLT,  SLATSQ , SLFS  and  TANDHI.  Also,  in  order 
1 to  provide  a basis  for  a necessary  ionospheric  model  we  have  chosen  to  use 

Chius  9 Fortran  computer  codes  IONDEN,  POLAR,  PSI,  SEMIAN,  TVARFZ,  TVEF1,  W, 

YONII  and  ZETA  to  compute  the  lower  ionospheric  density  (R*f.  10).  Some 
modification  to  the  old  SPECTER  code  INJECT  or  SINJCT  were  also  required. 

Description  of  the  3 new  main  programs  and  the  21  new  subprograms,  together  with  t 

} i 

f flow  diagrams,  are  given  below.  ' 

j 

a.  Program  AMBCAP  [ 

i ! 

Purpose: 

This  program  computes  the  magnetic  coordinates  corresponding  to  the  geographic 
| burst  point,  the  ambient  capacitance  and  the  height  integrated  ionospheric 

J resistance  along  the  field  line  passing  through  the  burst  point.  The  orthogonal 

j magnetic  grid  system  used  by  this  program  and  later  by  DEBCON  to  compute  the  j 

| debris  tube  motion  is  also  computed.  See  Figure  25. 

1 ! 


Labelling  information  for  each  case 
Default  option  index 

Altitude  of  effective  bomb  center 
Geographic  latitude  of  bomb  center 
Geographic  East  Longitude  of  bomb  center 
Local  time  in  hours  from  midnight 
Detonation  day  of  the  year 


Input : 

HCOM 

NUSE 

ALT 

DLAT 

DLONG 

TL 

DATE 
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f , NIS  No.  of  integration  steps 

| TEMP  Ionospheric  temperature 

1 1 BIGF  Instantaneous  10.7  cm  solar  flux 

> FBAR  10.7  cm  solar  flux  averaged  for  one  month 

H 

HM  Altitude  of  atmospheric  cutoff 

ELMAX  Maximum  L value  for  L grids 

RZUR  Smoothed  Zurick  sunspot  number 

i 

Output: 

Output  from  this  program  is  to  TAFE10. 

RECORD-1  RE,  BE,  SQ3,  ZMIN,  NIS,  NISMAX,  MODE,  DMLNG,  EL,  REQ,  OHM,  BDCK, 

CK(250,2),  PK(250,2).  See  the  AMBCAP  code  for  definition  of  variables). 

RECORD-2  F(100),  100  word  array  of  useful  Fortran  variables  (initial  file  A 
header  Information). 

Externals : 

ACATL,  ATMOS,  BLIMIT,  CONMAG,  DPLATR,  DTEQU,  EOF,  HEDGE,  INVAR,  MAGDEN,  SINMLT 
|{  and  VALUE. 

Method: 

Subroutine  INVAR  is  used  to  establish  the  burst  point  magnetic  coordinates  of  B, 

L and  the  magnetic  latitude.  Subroutine  MAGDEN  is  used  to  compute  the  geo- 
magnetic longitude.  The  ambient  capacitance  and  the  height-integrated  ionospheric 
resistance  are  computed  as  described  in  Section  II. 2b.  The  polar  intercept  array 
that  determine  the  orthogonal  magnetic  computational  grid  is  computed  by  speci- 
fying equal  altitude  intervals  between  the  atmospheric  cutoff  and  the  equator 
of  the  maximum  chosen  field  line  through  which  the  orthogonal  curves  must  pass 
(See  Figure  15). 


AMBCAP 


COMPUTE  J&C 
B*  & L*  VALUES 
AT  BURST  POINT 


WITH  AZ  a ZMAX 
COMPUTE  POLAR 
INTERCEPTS  AT 
EACH  ALTITUDE 
FOR  L = L 


COMPUTE  DIPOLE 
\*m  & VALUES 

AT  BURST  POINT 


PRINT 

RESULTS 


NORMAL' 

STOP 


Figure  25.  Program  AMBCAP 
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b.  Program  DEBCON 

\ 

Purpose : 

This  program  reads  in  the  AMBCAP  values  of  ambient  capacitance,  ionospheric  j 

resistance,  burst  magnetic  coordinates  and  the  polar  intercept  array  together 

with  necessary  card  input  values  specifying  the  debris  flux  in  order  to  compute  j 

the  motion  of  the  debris  tube  and  the  resulting  debris  density  distribution 

along  the  flux  tube.  The  debris  distribution  along  the  tube  is  determined  | 

from  the  motion  of  the  tube,  the  given  velocity  and  fission  flux  vs.  time  pro-  j 

i 

files  at  reference  planes  above  and  below  the  burst  point  and  the  fission  and  j 

debris  yield  of  the  burst.  See  Figure  26.  ' i 

i 

| 

Input : j 

TAPE10  See  output  of  program  AMBCAP  1 ] 

| 

HCOM  Labelling  information  for  each  case 

NUSE  Default  option  index  j] 

'i 

i j 
: i 

NT  No.  of  times  where  the  fission  debris  velocity  and  flux  are  specified 

YLD  Weapon  yield  in  megatons  ; j 

WF  Fraction  of  yield  that  is  fission  j 

RB  Radius  of  burst  bubble 

DSP  Initial  distance  of  upper  reference  frame  from  the  burst  point 

DSN  Initial  distance  of  lower  reference  frame  from  the  burst  point 

ETA  Time  independent  ratio  of  fission  fragment  to  debris  ions 

AMI  Time  independent  mass  per  debris  ion  ■ ■ 

TF  Maximum  time  the  motion  of  the  debris  tube  is  computed 

DT  Time  increment  used  by  the  code  to  compute  the  tube  motion 

R Ionospheric  resistance  (if  other  than  AMBCAP  value  is  to  be  used) 

TD  The  array  of  NT  times  where  the  fission  debris  distribution  is  specified 

VHP  The  array  of  NT  velocities  parallel  to  the  field  line  for  the  upper  ] 

reference  plane  corresponding  to  the  times  in  TD 
FP  The  array  of  fission  fragment  flux  at  the  upper  reference  plane  for 

the  NT  times  in  TD 
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VIIN  Same  as  VHP  but  for  lower  reference  plane 
FN  Same  as  FP  but  for  lower  reference  plane 

ETAT  The  array  of  the  ratio  of  fission  fragments  to  debris  ions  for  each 
of  the  NT  times  in  TD 

AMIT  The  mass  per  debris  ion  at  the  NT  times  in  TD 

NSBP  No.  of  computational  time  steps  between  print  outs  of  the  tube  position 

NDLG  No.  of  L values  to  output  to  file  A 

NTT  No.  of  times  per  L value  to  use  in  the  file  A output. 

NTPI  No.  of  DEBATL  print  control  times  to  read  into  array  TPI 

TPI  An  array  of  DEBATL  print  control  times  where  the  odd  values  specify 

time  print  increments  to  use  between  the  even  values  of  TPI 
(TPI(O)  = 0.  is  assumed) 

Output: 

Output  of  this  program  is  to  TAPE11  and  TAPE18 


TAPE11 
RECORD- 1 


RECORD-2 


The  no.  of  polar  intercepts 
The  width  of  the  flux  tube 
The  array  of  polar  intercepts 


1st.  L-shell  value 

Time  when  the  tube  is  at  this  L-shell 

Array  of  fission  debris  density  values  at  each  of  the 

NIS  points  and  each  hemisphere  for  this  L-shell 


RECORD-NDLG+1  EL 


Last  L-shell  value 

Time  when  the  tube  is  at  this  L-shell 

Array  of  fission  debris  density  values  at  each  of  the  NIS 
points  and  each  hemisphere  for  this  L-shell 


TAPE18  File  A (See  output  of  DEBRIS  code) 
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DEBCON 


READ  AND 
PRINT  CARD 
INPUT 


INPUT 


VALID 


INITIALIZE 
VALUES  . 
T,  L,  Q,  Q 


T = TMAX 

X 

L > L 


COMPUTE 
TUBE  MOTION 
AT  TIME 
T = T + DT 


COMPUTE 
F.F.  DISTRIBUTION 
AT  TIME 
T = T + DT 


L < L 


COMPUTE 
FILE  A 
DATA 


COMPUTE  Q, 
THE  INITIAL  ' 
TOTAL 
CURRENT 


'OUTPUT^ 

F.F. 

S.DATA  a 


f OUTPUT^ 

( NORMAL\ 

V FILE  a 1 

ySTOP  I 

Figure  26.  Program  DEBCON. 
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Externals: 

ACATL,  ARANGD , DEBATL,  DRIFTC,  HEDGE,  MOVEL,  SINDEX  and  VALUE 


!i 


Method: 

This  program  computes  the  position  of  the  flux  tube  and  the  fission-fragment 
distribution  profile  within  the  flux  tube  at  discrete  times.  The  tube  motion 
is  governed  indirectly  by  the  equivalent  RLC  circuit  equation  (Eq.  43). 

The  height  integrated  ionospheric  resistance  is  provided  by  AMBCAP  or  is 
defined  via  card  input.  Its  value  is  assumed  to  remain  constant. 

At  each  time  step  an  inline  function  routine  FLF  computes  the  inductance  from 
Eq.  (17),  the  ambient  capacitance  is  computed  by  subroutine  ACATL  while  the 
total  capacitance  and  azimuthal  drift  current  are  computed  by  subroutine  DBBATL 
which  also  updates  the  fission  fragment  distribution  profile.  The  azimuthal 
drift  current  and  RLC  circuit  values  are  used  by  subroutine  MOVEL  to  solve  the 
circuit  equation  and  thus  obtain  the  net  charge  across  the  tube.  Thus, 

^ can  then  be  computed  from  Eq.  (42)  and  the  new  position  of  the  flux  tube 
evaluated.  The  above  steps  are  repeated  until  the  flux  tube  exceeds  given 
L-shell  bounds  or  until  the  time  exceeds  the  maximum  allowed. 


c.  Program  DSTRIB 

Purpose: 

This  program  computes  the  redistribution  of  injected  electrons  within  the  flux 
tube  due  to  the  generated  electric  field  that  moves  the  flux  tube.  The  re- 
distribution is  computed  for  electrons  at  the  eastern  edge  of  the  flux  tube. 

The  energy  distribution  of  the  redistributed  electrons  at  the  equator  is  also 
computed.  This  is  necessary  since  these  electrons  will  in  general  no  longer 
have  the  simple  fission  beta  energy  spectrum.  This  complication  in  the  spectrum 
is  due  to  the  fact  that  electrons  with  different  energy  do  not  drift  to  the  east 
at  the  same  rates. 
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Input : 

TAPE 11 

DTMAX 
DU 
FK 

NEM 

EG 

T 

PHIED 

Output: 

At  present  there  is  no  output  tape  from  this  program.  However,  the  printed 
output  provides  details  of  the  energy  and  pitch  angle  distribution  of  the 
redistributed  electrons. 

Externals: 

BLIMIT,  BOUNCE,  BSFRMP , EXPINT  LESSEN,  LINTRP,  PHIDOT,  SINDEX,  SLATSQ , VALUE 
and  YINOX. 

j 

' Method : 

The  method  used  by  this  program  is  described  in  Section  n.3a.  The  use  of 
Liouville's  equation  together  with  conservation  of  the  first  and  second  adiabatic 
invariants, conservation  of  magnetic  moment,  and  a uniform  source  in  longitude 
result  in  Eq.  108  that  relates  the  time  derivative  of  the  directional  flux 
^ to  a double  integral  over  the  fission  fragment  source  function  along  the 
appropriate  field  line  at  each  possible  longitude  within  the  flux  tube. 

This  flux  ^ gives  the  number  of  electrons  per  square  cm  per  sec.  per  sec. 
per  sterladian  per  MeV  at  the  eastern  edge  of  the  flux  tube.  This  equantity 
can  be  Integrated  over  time  within  the  flux  tube  and  over  pitch  angle  to  obtain 


See  output  of  program  DEBCON 
Maximum  time  step  increment  allowed 

The  cosine  (pitch  angle)  increment  to  use  for  the  pitch  angle  grid. 
The  number  of  longitude  intervals  to  use  within  the  flux  tube. 

The  number  of  discrete  energy  values  where  the  electron  distribution 
is  to  be  computed. 

The  array  of  NEM  energy  values 

The  time  when  the  redistribution  is  to  be  computed  (the  time  T 
specified  the  L-shell  position  of  the  flux  tube) 

The  east  longitude  of  the  eastern  edge  of  the  debris  tube. 
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the  desired  energy  distribution  at  the  equator.  A further  integration  over 
energy  results  in  an  omnidirectional  flux  that  when  divided  by  2it/A4> 
can  be  compared  with  the  equatorial  longitude  independent  flux  that  is  computed 
by  other  SPECTER  codes. 


d.  Subroutine  ACATL 


Purpose: 

This  subroutine  computes  for  a given  L-shell  the  ambient  capacitance  of  the  flux 
tube  (assuming  the  frozen-field  model)  the  distance  to  the  magnetic  equator 
from  the  earth's  center,  the  distance  to  the  equator  along  the  field  line  from 
the  point  of  intersection  of  the  orthogonal  curve  defined  by  each  polar  intercept 
and  the  magnetic  field  intensity  corresponding  to  these  intersections. 


Input : 
ARG 


ACATLC 


Output: 

AC 

REQ 


SK 


BK 


FLATK 


Task  option,  set  to  zero  if  ambient  capacitance  and  distance  to  equator 
REQ  are  not  required 


Labeled  COMMON  variables  RE,  ZMIN,  KM,  MODE,  DMLONG,  EL,  BDCK, 
CK  and  PK  (See  ACATL  code  for  definitions  of  variables) 


Ambient  capacitance  of  the  flux  tube  at  the  L-shell  EL  in  farads 
Distance  in  earth  radii  from  the  center  of  the  earth  to  the  magnetic 
equator 

An  array  of  distances  in  KM  along  the  field  line  EL  to  the  equator 
from  the  intersection  points  defined  by  the  polar  intercepts  in  PK 
for  each  hemisphere. 

An  array  of  magnetic  field  intensities  in  gauss  along  the  field  line 
EL  at  the  intersection  points  defined  by  the  polar  intercepts  in  PK 
for  each  hemisphere 


2 h 

An  array  of  magnetic  latitude  parameters  equal  to  (1+3  sin  A ) * 


where  A is  of  the  magnetic  latitude  at  each  intersection  point  defined 
by  the  values  of  polar  intercepts  in  PK  for  each  hemisphere. 
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Subroutine  BLIMIT  is  called  to  define  the  equatorial  value  of  magnetic  intensity 

B at  the  input  L-shell.  For  each  polar  intercept,  subroutine  DPLATR  is  then 
o 

called  to  compute  the  magnetic  latitude  that  defines  the  magnetic  field 
intensity  since 

B = Bq  (1  + 3 sin2  \ )*5/cos6X  . (123) 
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The  function  subroutine  DTEQU  is  called  to  compute, as  a function  of  the  distance 
along  the  field  line,  a quantity  used  later  by  BEBATL  to  determine  the  debris 
distribution.  Subroutined  CONMAG  is  called  with  the  magnetic  latitude  and  input 
longitude  as  arguments  in  order  to  compute  the  geographic  altitude  of  the  inter- 
section points  that  define  the  flux  tube  segments.  The  contribution  to  the 
ambient  capacitance  from  each  flux  tube  segment  k is  then  computed  from  the 
equation 


"ko 


<\  + + 


(124) 


where 


Cko  • »ko  aSko/Bko'  \ ■ « + ’ si”-»  k^ 


and  the  o subscript  refers  to  the  initial  position  of  the  flux  tube  where  the 
mean  ionospheric  density  in  the  segment  was  p^0  and  the  segment  length  was  AS^. 
The  above  equation  is  used  provided  the  altitude  of  the  segment  is  above  the 
atmospheric  cutoff.  An  interpolation  is  made  if  the  tube  segment  lies  partly 
above  and  below  the  atmospheric  cutoff.  The  total  ambient  capacitance  is  obtained 
by  summing  over  each  C^.  The  distance  to  the  equator  from  the  earth’s  center 
is  computed  from  the  last  altitude  computed  (corresponding  to  X=o). 
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e.  Subroutine  ARANGD 


Purpose: 

Subroutine  ARANGD  reads  in  and  arranges  a subset  of  the  computed  debris 
distribution  data  on  TAPE11  into  the  format  specified  by  file  A and  writes  the 
arranged  data  on  TAPE18.  This  data  can  then  be  processed  by  the  SEPCTER 
injection  and  flux  codes.  Only  a subset  of  the  DEBATL  generated  data  is  written 
to  TAPE11,  however,  this  data  in  general  is  specified  for  many  more  points  along 
a field  line  than  is  practical  for  output  as  file  A and  must  be  reduced.  Also, 
since  the  tube  motion  may  be  sinusoidal  in  L,  data  on  TAPE11  from  various  time 
segments  must  in  general  be  combined  in  order  to  specify  the  debris  at  any  given 
L-shell  as  a function  of  time  as  is  required  by  the  file  A format. 


Input : 
NLG 

ELG 

NLT 

ELTBL 

TI 

TF 

DT 

CDEL 

NT 

ITL 

ITD 


The  number  of  L values  in  array  ELG  where  the  fission  debris  density 
is  to  be  evaluated 

The  array  of  monotonic  increasing  L-shell  values  that  covers  the  range 
of  tube  motion.  The  debris  is  to  be  spedified  at  these  L values 
The  number  of  L values  in  array  ELTBL  and  the  number  of  field  line 

debris  distribution  density  profiles  written  to  TAPE11 

< 

The  array  of  L values  where  the  debris  distribution  density  profiles 
are  contained  on  TAPE11.  These  L values  are  separated  by  a constant 
time  increment. 

The  initial  time  relative  to  the  burst  when  data  is  first  written  to 


TAPE 11  (TI  - 0). 


The  last  time  when  data  is  written 

The  time  increment  in  seconds  between  the  stored  data  on  TAPE11 

The  ratio  at  the  equator  of  the  height  of  the  debris  tube  to  the 

distance  to  the  equator  from  the  earth's  center 

The  number  of  times  when  the  field  line  debris  density  profiles 

are  to  be  computed  for  each  L value 

The  Input  tape  number  (11) 

The  output  tape  number  for  file  A (18) 
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Output : 


ARANDG  writes  data  on  TAPE18  according  to  the  format  described  below: 

Record  No.  Description 

2 A monotonic  increasing  array  of  L values  for  which  the  density  is 
specified. 

3 An  array  that  specifies  the  number  of  points  used  along  each  L shell. 

4 An  array  specifying  distance  to  the  equator  in  km  for  each  point. 

5 An  array  specifying  magnetic  intensity  in  gauss  for  each  point. 

6 An  array  specifying  the  corresponding  equatorial  cosine  (pitch 

andle)  for  each  point. 

7 An  array  specifying  the  time  for  each  point  at  which  the  debris 
number  density  is  given  in  the  following  record. 

8 An  array  specifying  the  debris  number  density  for  each  point 
corresponding  to  the  times  specified  in  the  previous  record. 


The  records  between  the  dashed  lines  are  repeated  for  every  L value  contained 
in  record  2 and  the  records  between  the  dotted  lines  are  repeated  for  as  many 
times  as  is  adequate  to  represent  the  desired  debris  distribution. 

Method: 

At  the  present  time  ARANGD  handles  as  many  as  six  entries  of  the  debris  tube 

into  a particular  L interval.  Entrance  and  exit  times,  t^  and  t , are 

in  out 

computed  for  the  times  when  the  center  lines  of  the  debris  tube  cross  the  L values, 

L . and  L , which  are  separated  approximately  by  the  tube  width  L.  The 
min  max 

L limits  are  given  by  the  equations. 


L . + 

min 


where  L is  the  location  of  the  tube.  Since  the  tube  width,  L,  is  proportional 
2 

to  L , the  above  equations  can  be  solved  to  give. 
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(127) 


Lmin  " 1-1  + + 2AL/L)1*]  L2/  L 

Lmax  " [1_  (1  “ 2 ^L/L)1*]  L2/  L 

i L * (L/L  )2  A L 
o o 

and  A L * L Adi 
o o 

The  debris  density  for  the  specified  L values  can  then  be  computed  by  inter- 
polation at  the  times  t i t £ t *ro®  the  previously-computed  debris  distri- 
bution data  on  input  TAPE11. 


(128) 

(129) 
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f.  Function  BOUNCE 

Purpose; 

This  routine  computes  the  electron  bounce  period  as  a function  of  energy,  pitch 
angle  and  L-shell. 

Input : 

EL  The  L-shell 

EMEV  The  energy  in  meV. 

S The  value  1 sine  (pitch  angle) 

$ 
l 

Output:  | 

BOUNCE  The  electron  bounce  period  in  seconds  j 

] 

Method:  I 

The  bounce  period  T is  computed  from  the  equation  I 


4 


i 


»>.(!•:)  - (l  + K/n/|(i  + h/ iO 2 - 1 l\  (l'i'i) 

h(S)  - 1.180171  - ().61%9IS0,7\  (134; 

K = .51098,  L = EL,  K = KMKV  and  T = BOUNCE 

i> 

Tho  routine  attempts  to  economize  by  not  recomputing  f (L) , g(E)  or  h(S)  if 
the  corresponding  input  values  do  not  change. 


g.  Subroutine  BSFRMP 


« 
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Purpose: 

This  routine  computes  for  a given  field  line  an  array  of  magnetic  intensities 
and  distances  along  the  field  line  to  the  equator  corresponding  to  the  inter- 
section points  defined  by  an  input  array  of  polar  intercepts. 

Input : 

EL  The  magnetic  L-shell  or  field  line 

KM  The  number  of  polar  intercepts  or  intersection  points 

KMAX  The  1st.  dimension  of  arrays  PK,  BK  and  SK 

PK  The  array  of  polar  intercepts  in  earth  radii  for  each  hemisphere 

Output: 

BK  The  array  of  magnetic  intensities  in  gauss  for  each  hemisphere 

SK  The  array  of  distances  along  the  field  line  to  the  equator  for  each 

hemisphere.  Distances  are  negative  in  the  southern  hemisphere. 


Externals: 

BLIMIT,  DPLATR  and  DTEQU 

Method: 

The  equatorial  magnetic  intensity  Bq  at  L-shell  L is  computed  by  calling  subroutine 
BLIMIT.  For  each  point  on  the  field  line  defined  as  the  intersection  of  the  ortho- 
gonal curve  that  has  an  altitude  equal  to  the  polar  intercept  values  at  the  poles, 
subroutine  SPLATR  is  called  to  compute  the  magnetic  latitude  . The  magnetic  in- 
tensity B is  computed  from  Eq.  123.  Subroutine  DTEQU  is  then  called  with  L and 
as  input  arguments  to  compute  the  distance  to  the  magnetic  equator. 
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h.  Subroutine  CONMAG 


/ 


ft 


S-V'. 


Purpose : 

Subroutine  CONMAG  is  used  to  convert  from  geographic  coordinates  to  magnetic 
coordinates  when  the  accuracy  required  (between  the  real  magnetic  field  and 
geographic  position)  is  not  larger  and/or  where  the  computer  time  for  such 
conversions  would  be  prohibitive  otherwise  (a  modified  version  of  INVAR 
could  be  used  for  this  purpose  if  great  accuracy  was  required  and  computer 
time  of  no  concern).  This  routine  is  needed  since,  in  general,  the  neutral 
and  ionospheric  model  atmospheres  required  in  this  code  use  geographic 
coordinates  as  independent  variables  whereas  magnetic  coordinates  must  be 
used  to  describe  charged  particle  motion. 


Input ; 

M A mode  option  index  (M=l,2,3,4,5,  or  6) 


M=l,3  or  5 for  centered  dipole  approximation 
>^2,4  or  6 for  displaced  dipole  approximation 


ALT 

Geographic  altitude 

M-l ,2,3,  or  4 

DLAT 

Geographic  latitude  in  degrees 

M-l  or  2 

DLONG 

Geographic  longitude  in  degrees 

M=1  or  2 

DMLAT 

Sign  of  geomagnetic  latitude 

M=3  or  4 

DMLONG 

Geomagnetic  latitude  in  degrees 

Geomagnetic  longitude  in  degrees 

M-5  or  6 

M-3 ,4,5  or 

6 

EL 

Magnetic  L-shell  value 

M-3, 4, 5 or 

6 

Output : 

ALT 

Geographic  altitude 

M-5  or  6 

DLAT 

Geographic  latitude  in  degrees 

M-3  or  4 

DLONG 

Geographic  longitude  in  degrees 

M-3  or  4 

DMLAT 

Geomagnetic  latitude  in  degrees 

M-1,2,3  or 

4 

DMLONG 

Geomagnetic  longitude  in  degrees 

M-l  or  2 

DDIP 

Magnetic  dip  angle  in  degrees 

M-1,2,3  or 

4 

B 

Magnetic  intensity  in  gauss 

M-1,2,3  or 

4 

EL 

Magnetic  L-shell  value 

M-l  or  2 
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Method: 


For  modes  M«1  or  2 geographic  cartesian  coordinates  Xg  ■ (Xg.Yg.Zg)  are 
computed  from  the  equations 


= (1  + h/R£)  cos  cos  $g 

(135) 

= (1  + h/Rg)  cos  XG  sin  <J>G 

(136) 

* (1  + h/RE)  sin  XG 

(137) 

where  h is  the  altitude,  R^  is  the  earth's  radius,  XG  is  geographic 
latitude  and  $G  is  geographic  longitude.  These  coordinates  are  then 
translated  to  a magnetic  origin  by  the  matrix  equation 


*T  " XG  + T 


(138) 


where  T = (.05764855,  - .03208677,  - .01842101)  if  M = 2 and  T = 0 

if  M * 1.  The  magnetic  cartesian  coordinates  X are  then  computed  by 

m 

the  rotation 


X 

m 


= AX 


T 


where 


A 


( 


.97129205 
.23641749 
- .02642934 


.22646696 

.95294113 

.20153391 


.07283175 

.18976292 

.97912490 


(139) 


(140) 


and 


The  desired  spherical  magnetic  coordinates  are  then  computed  using  the 
relations 


A 

m 


sin"1  (Z/R  ) 
m m 


(141) 


4 


•„  — JU- 
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(142) 


For  modes  M * 3 or  5,  the  X coordinates  are  first  computed  and  the  above 

m 

equations  used  in  an  iterative  way  to  compute  the  corresponding  geographic 

position  (we  use  the  fact  that  the  dipole  displacement  T is  small  compared 

with  the  distances  to  regions  outside  the  earth  to  compute  R^) . For  modes 

M = 5 or  6,  the  X coordinate  is  also  computed  first  given  R . The  geograph- 
m in 

ic  coordinates  are  then  computed  as  above  from 


I 


L 


X„  - A_1x  - T. 

Vj  TU 


The  altitude  can  then  be  obtained  since 


(146) 


r 2 2 2 "1*^ 

h = 6371.2  |^(XG  + Yg  + ZG  ) - 1J  . 


(147) 


i.  Subroutine  DEBATL 

Purpose : 

Subroutine  is  used  to  perform  three  basic  tasks.  First,  it  is  used  to 
compute  some  required  burst  parameters  and  to  print  out  the  input  debris 
distribution  data  for  each  data  plane.  Secondly,  it  is  used  to  compute 
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polar  intercepts  for  the  burst  point  and  the  data  planes,  to  establish 
indicies  relating  the  position  of  the  data  planes  in  the  computational 
grid,  to  compute  the  total  number  of  fission  fragments,  the  total  hot- 
ion  mass  and  energy  and  then  to  scale  the  input  distribution  to  the  total 
fission  yield  as  given  by  the  input  data  (these  quantities  should  be  con- 
sistent but  may  not  be!).  Finally,  DEBATL  is  used  to  compute  the  fission 
fragment  density  distribution  profile  along  the  length  of  the  flux  tube 
as  is  dictated  by  the  motion  of  the  tube  on  the  flux  and  velocity  distribu- 
tion at  the  data  planes. 


Input : 


Time  in  seconds 


ACATLC  The  labeled  common  array  variables  RE,  BE,  NIS,  DMLONG,  EL,  BDCK, 
PK,  SK,  BK  and  FLATR 

DEBC  The  labeled  common  array  variables  TI,  ELI,  IH,  BX,  DPR,  C,  CLDG, 
Q,  BALT,  BPHID,  BLATD,  YLD,  WF,  BR  and  DEBIN 

DATAC  The  labeled  common  array  variables  NT,  DSP,  DSN,  NTPI,  TPI,  ETAT, 
AMIT,  TD,  VHP,  FP,  VIIN,  FN  and  DATA 
(See  DEBATL  code  for  definitions  of  variables.) 


Output : 

TAPE9  Print  file  (should  be  rewound  and  copies  to  OUTPUT  to  be  listed) 

DEBC  The  labeled  common  array  variables  DPHI,  SET,  SE,  BEQ,  BY,  C, 

ELO,  VI,  TDEB , TTBP  and  SIGMA 

DATAC  The  labeled  common  array  variables  FP  and  FN 

MOVEC  The  labeled  common  array  variables  ADI  and  ADID 

DEBDAT  The  labeled  common  array  variable  EDEN 

(See  DEBATL  code  for  definitions  of  variables.) 


Externals; 


BLIMIT,  DPLATR,  DTEQU,  LINTRP,  SFROMB,  SLFS,  XINDX  and  YINDX. 
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Method: 

DEBATL  defines  the  burst  area  A_  in  terms  of  the  bubble  radius  r,  where 

o D 


Afi  = 711 


The  constant  azimuthal  tube  width  is  defined  to  be 


(148) 


= /V(RELb  cos  V 


(149) 


where  is  the  radius  of  the  earth,  is  the  burst  L-shell,  and  X^  is  the 
magnetic  latitude-  The  total  number  of  available  fission  fragments,  n^t  is 
computed  assuming 


nft  = 1.  x 10  Wf*YLD 


(150) 


where  Wf  is  the  fission  fraction  and  YLD  is  the  yield  in  megatons. 

However,  the  total  number  of  fission  fragments  Nft  is  also  defined  by  the 
fission  fragment  flux  at  the  input  data  planes  since 


A-  A. 

[A+jF+  dt  + A_  jF_  dt] 


(151) 


where  A. 


is  the  initial  area  of  the  higher  data  plane  where  the 


magnetic  intensity  is  B+,  A+  = A^B^/B_  is  the  initial  area  of  the  lower  data 
plane  where  the  magnetic  intensity  is  B_,  F+  is  the  fission  fragment  flux 
at  the  + data  plane,  F_  is  the  fission  fragment  flux  at  the  - data  plane,  T 
is  the  maximum  time,  and  B^  is  the  magnetic  intensity  at  the  burst  point.  Thus, 
this  code  defines  a scaling  coefficient  Cff  such  that 


Cff  “ nff/Nff 


(152) 


and  new  flux  values  f+  = F+  and  f_  *■  F_  so  that  the  total  number  of 
fission  fragments  is  leaving  the  data  planes  is  n^f. 
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The  flux  tube  is  broken  up  into  three  distinct  regions  (region  1,  2 and  3) 

r 

for  the  purpose  of  computing  the  fission  fragment  profile  along  the  flux 
, tube.  Region  1 is  defined  as  the  portion  of  the  flux  tube  below  the  lower 

data  plane  (this  region  exists  in  only  one  hemisphere);  region  2 is  defined 
as  the  portion  of  the  flux  tube  between  the  lower  data  plane  and  the  upper 
data  plane  (this  region  includes  the  burst  point).  Region  3 is  defined  as 
the  remaining  portion  of  the  flux  tube  (this  region  can  exist  in  both  hemis- 
j pheres) . The  fission-fragment  distribution  in  region  2 is  ambiguous  and  we 

have  chosen  to  compute  the  unambiguous  mean  density  of  fission  fragments  in 
I this  region  and  assume  it  is  valid  everywhere  within  the  region.  Regions 

1 and  3 are  similar  except  that  region  1 is  constrained  to  the  burst  hemis- 
phere . 
i 

The  fission  fragment  number  density  at  time  t in  region  2 is  computed  from 
the  equation 


n2(t)  = n2(t)  " K(t)  + ft(t,)  + f^(t)  + At' )]--(■-)  (153) 

2V2 

i 

I where  t'  = t - At  and  V„  is  the  volume  of  the  flux  tube  in  region  2, 

i a A z A A 

j f+  * A+f+  and  f_.  = A_f_.  (Note  that  f and  f_  are  the  number  of  fission 

fragments/sec . Leaving  the  data  planes  and  that  these  are  conserved 
quantities  independent  of  the  motion  of  the  tube.) 

i 

In  regions  1 and  3,  the  debris  distribution  is  assumed  to  have  a front 
■:  corresponding  to  time  0 when  particles  first  left  the  data  plane  and  a 

tail  corresponding  to  time  t if  t<T  or  time  T if  time  t>T  when  particles 
J last  left  the  data  plane.  At  each  time  step  and  for  each  point  along  the 

flux  tube  occupied  by  debris,  an  array  of  times  t and  velocities  V are 
| kept  that  indicate  when  the  particles  at  the  point  in  question  left  the 

[ data  planes  and  the  velocity  they  had  when  they  left.  Similar  times  and 

velocities  are  kept  for  the  head  and  tail  of  the  distribution.  Thus,  at 
each  time  step  the  new  position  of  the  head  and  tail,  as  well  as  that  of 
the  particles  that  were  at  the  given  grid  points,  can  be  computed  from  the 


90 


old  velocity,  V.  The  arrays  t and  V can  then  be  updated  by  interpolating 
in  distance  along  the  field  line  (note  that  we  assume  that  the  velocity 
does  not  change  in  the  particles  frame  of  reference).  Then,  for  each 
grid  point  k between  the  front  and  tail  f_^(rk)  say,  can  be  computed 
from  which  the  number  density  can  be  expressed  as 


pff(k,t)  * fJ(Tk)(Tk  - Tk]/|A+(k)*(Sk  - Sk)] 


(15«) 


where  rk  - tk(t  - At),  = Sk(t  - At)  and  A+(k)  is  the  cross-sectional 
area  of  the  flux  tube  at  point  k and  time  t. 


j • Subroutine  DPLATR 


Purpose : 


This  subroutine  computes  the  dipole  magnetic  latitude  X at  the  point  where 

m 


the  L-shell  defined  by 


R = L cos  X 


intersects  the  orthogonal  shell  defined  by  the  relation 


(155) 


R * P/sin1/2X 


(156) 


Input : 


The  magnetic  L-shell 
The  polar  intercept 


Output: 


RLAT  The  magnetic  latitude  X in  radians. 

m 


Method : 


The  subroutine  DPLATR  uses  the  following  method  to  compute  Xk  given  P^  and  L. 
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is  computed  using  the  equations: 


U 


T 


First,  an  approximate  value,  X°  , 


,(u./(2.9  u,  + 1.5)) 


V1  • -41 V 

(L°810  Uk 
n/2  - [ ( 1— . 25/  /uk)/  /uk 


% 


, Pk  * 0 
, uk<  1.3 

> uk  75»  \i  1-3  (157) 


2 

where  u^  * (L/pk)  . 

The  approximate  value,  X°,  is  then  refined  by  using  Newton's  method 
with: 

Xk+1  = Xk  ' fi/(df/dXk)i»  i = 

where  = f^X^)  = uk  cos4  X^  - sin  X^  (158) 

and  \ = -cos  X^  (4u^  cos^  X^  sin  X^  + 1).  (159) 

' k'j.  & 

The  iteration  sequence  is  stopped  when  | (A^+1)Ak  | < e * 10  6. 


k.  Function  DTEQU 

Purpose : 

This  routine  computes  the  distance  along  a dipole  field  line  to  the  equator. 


Input : 

EL  The  L-shell  defining  field  line 

X The  quantity  sine  (magnetic  latitude) 

1/2 

Y The  quantity  (1+3x2) 

Output: 

DTEQU  The  distance  to  the  equator  in  km  where  distances  are  positive 
north  to  south. 
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Method : 

The  distance  to  the  equator  DTEQU  is  evaluated  using  the  equation 

DTEQU  = sign  [3185.6  EL(|x.|y  + ln(/T  |x|  + y)//3),  xj  (160) 


1.  Subroutine  EXPINT 


This  subroutine  computes  the  integral 
Xo 


i-Jr 


(161) 


Input : 


XI  The  value  of  the  independent  variable  at  the  lower  limit  of 

integration 

Y1  The  value  of  the  integrand  at  XI 

X2  The  value  of  the  independent  variable  at  the  upper  limit  of 

integration 

Y2  The  value  of  the  integrand  at  X2 

Output: 

DINT  The  positive  value  J I [ 

A The  reciprocal  e folding  distance  between  XI  and  X2 


Method : 

Y1>0 

DINT  = (Y2  - Yl)/A;  if  Y2>0 

Y1£Y2 


(162) 


where  A = In  (Y2/Yl)/(X2  - XI), 


(163) 


and  DINT  - (Y2  + Yl) j (X2  - Xl)|/2 


(164) 


otherwise 
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m.  Subroutine  IONCON 


Purpose : 

This  routine  interpolates  in  altitude  from  a table  of  observed  relative 
ion  concentrations  vs.  altitude  below  500  km,  to  compute  the  ion  concen- 
trations of  0+,  NO"*",  0+,  Hg+  and  H+  given  the  total  ion  concentration 
and  altitude. 


Input : 


The  altitude  where  the  ion  concentration  is  TC 


TC  The  total  ion  concentration  or  number  density  at  altitude  H 


Output : 


Cl  An  array  of  0^,  N0+,  0+,  He+  and  H+  ion  concentration  at  altitude  H. 


Externals: 


XINDX 


Method : 

The  relative  ion  concentrations  are  interpolated  from  a table  taken 
from  the  results  of  C.  Y.  Johnson  (Ref.  11).  The  above  ion  number  densities 
n^  are  then  computed  from  the  equation 


n±  =*  (N/I  Ri>  R±, 


(165) 


where  N m TC,  n^  * n(02(»  = n(N0+),  . . . , n^  = n(H+)  and  n^  = CI(i). 


n.  Subroutine  LESSEN 


Purpose : 


This  subroutine  is  used  to  store  a subset  of  points  from  one  array  into 
another  array  each  of  which  are  two  dimensional  with  the  second  dimension 
equal  to  two  corresponding  to  two  hemispheres. 
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Input : 

NI  The  number  of  val ues/hemisphere  in  AT 

NID  The  first  dimension  of  AI 

NFM  The  maximum  number  of  values/hemisphere  in  AI  to  be  placed  in  AF 
NFD  The  first  dimension  of  AF 

AI  The  array  of  values  to  be  chosen 


Output : 

NF  The  number  of  values /hemisphere  in  AI  placed  in  AF 

AF  The  array  of  placed  values  (note  that  AI  = AF  is  OK) 

Method : 

For  each  hemisphere  the  array  AF  is  formed  from  AT  such  that: 

AF(J)  = AI (K) , J-1,2 ,NF 

K=i . 1+2 i+(':f-i>  i 


(167) 


(168) 


where  I - 1 + (NI  - 1)/(NFM  -1) 


NF  = 1 + (NI  -1)/I 


o.  Subroutine  LINTRP 

Purpose : 

This  routine  does  linear  interpolation  over  a monotonic  increasing  array  of 
independent  variables. 

Input; 

M =>1,  when  independent  variable  has  changed 

*0,  when  X,  NX  and  XT  have  not  changed 
NX  The  number  values  in  arrays  XT  and  YT 

X The  value  where  Y is  to  be  elevated 

XT  The  independent  monotonic  increasing  variable  array 

YT  The  dependent  variable  array 
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The  interpolated  dependent  variable  corresponding  to  X 


I 

rj 


Output : 
Y 


Externals : 

XINDX 

Method: 

Straight  line  linear  interpolation  is  used.  Some  economy  in  computation 
can  be  obtained  as  in  DEBATL  where  several  quantities  are  interpolated  for 
the  same  value  of  independent  variable  by  setting  M=0,  in  which  case  the 
interpolating  coefficient  is  not  updated. 


p.  Subroutine  MAGDEN 


Purpose : 

A model  of  the  earth's  ambient  ion-plasma  is  required  to  evaluate  coeffi- 
cients in  the  differential  equation  that  governs  the  motion  of  the  debris 
tube.  MAGDEN  provides  a global  model  of  the  earth's  ambient  ion  concentra- 
tion as  a function  of  position  with  annual,  diurnal  and  solar-activity 
cycles. 


*1 


i 


Input : 


M =1,  for  centered  dipole  field 

*2,  for  displaced  dipole  field 

I 

ZL  Altitude  limit  of  CHIU's  IONDEN  code  in  km 

Z Altitude  of  geographic  point  in  km 

DGLT  Geographic  latitude  of  point  in  deg. 

DGLNG  Geographic  longitude  of  point  in  deg . 

RZUR  Zurich  smoothed  sunspot  number 

TL  The  local  time  in  hours 

TM  Time  in  months  from  December  15th  of  the  previous  year 
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dilt'iH  iighMrf"'-  


Output : 

C The  total  Ion  concentration  at  altitude  Z * 

Cl  The  array  of  02+,  N0+,  0+,  He+  and  H+  ion  concentrations  at 

altitude  Z 

DMLT  The  magnetic  latitude  in  deg. 

DMLNG  The  magnetic  longitude  in  deg. 

DDIP  The  dip  angle  of  the  field  line  in  deg. 

B The  magnetic  intensity  in  gauss 

EL  The  L-shell  value 

ELP  The  L-shell  value  of  the  plasmapause  at  time  TL 
TIK  The  ion  temperature  in  deg.  K. 

ZIKM  The  0+/H+  ion  transition  height  in  km  a 

j 

Externals:  ] 

■I 

CONMAG,  IONCON , IONDEN.  PPLOC  and  TANDHI  1 

Method:  | 



MAGDEN  computes  Ion  concentrations  by  extrapolating  as  required  Chiu’s 
ionosphere  model  (Ref. 10)  in  a theoretical  manner  that  is  consistent 
< with  observations.  The  ion  concentrations  are  extrapolated  beyond  alti- 

I tudes  of  500  km  with  the  help  of  new  subroutines  PPLOC,  TANDHI,  and  CONMAG. 

PPLOC  interpolates  to  compute  the  mean  L value  position  of  the  plasmapause 

as  a function  of  local  time.  TANDHI  computes  ion  temperatures  and  0+/H+ 

ion  transition  heights  as  a function  of  magnetic  latitude  and  local  time 

from  a model  based  on  the  work  of  Titheridge  (Ref.13).  CONMAG  computes 

the  magnetic  latitude,  Xm,  magnetic  longitude,  cc^,  field  intensity,  B,  and 

L- value  as  a function  of  geographic  position.  It  also  computes  geographic  j 

latitude,  magnetic  latitude,  and  intensity  for  a given  altitude,  field 

line  and  magnetic  longitude. 


i 


i 


l 

i • 

• I 
; ! 


IHWIIP  W- 


1.  The  L-value  of  the  plasmapaase  Lpp  is  computed  from  the  local  time. 

2.  The  magnetic  coordinates  (x^,  # )t  B,  and  L are  computed. 

3.  The  magnetic  latitude,  Xm  (1000)  at  an  altitude  of  1000  km  along  the 
field  line  (L,  com)  is  computed. 

4.  The  ion  temperature,  T,  and  0+/H+  transition  height,  h. , are  computed 
from  Xm  (1000),  local  time,  and  perhaps  season. 

5.  The  latitudes,  \m  (500)  and  \ (500),  at  an  altitude  of  500  km  along 

the  field  line  are  .imputed. 

6.  The  total  ion  concentration,  N_  - N (500),  is  computed  from  Chiu's 
model  where  we  assume  (500)  = fO  (500)]  + [He+(500)]  + [H  (500)]. 


i 7.  Assuming  diffusive  equilibrium  and  a constant  ion  temperature,  T,  the 

■ scale  heights  of  H+  and  0+  are  computed.  Then  the  ratio  [H+]/[0+]  at 

!•  i 

■ 500  km  is  determined  knowing  that  at  the  transition  altitude,  [0  (h^.)]  = 

j [H+(ht)1. 

8.  The  ratio  [H+]/[0+]  at  500  km  is  computed  by  using  the 
observation  that  at  this  altitude  [H+]/[He+]  « 7;  hence 
[He+]/ro+]  = [H+]/7ro+]. 

|{  9.  The  concentration  f0+(500)]  is  then  computed  from  the  equation 

j [0+]  = Nj/(l  + fH+]/[0+l  + [He+l/[0+]  whence  the  concentrations 

jl  [He+( 500)1  and  [H+(500)]  are  determined. 

For  50O  < h s ht  or  L 5 Lpp  

10.  The  concentrations  of  0+,  He+,  and  H+  are  computed  using  the  known 
values  at  500  km  (from  step  9)  and  the  computed  scale  heights. 

For  h > h^  and  L > Lpp 
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q.  Subroutine  MOVEL 

Purpose : 

This  subroutine  is  used  to  solve  the  equivalent  circuit  differential  Equation 
(43)  at  each  time  step  and  thus,  to  compute  the  motion  of  the  convecting  flux 
tube  via  Equation  (42). 


Input : 


Input  time  in  seconds 
Time  increment  of  step 
L-shell  value  at  input  time 

Charge  on  flux  tube  in  coulombs  at  input  time 
1st  time  derivative  of  Q at  input  time 
Ionospheric  resistance  in  ohms 

Capacitance  of  flux  tube  in  farads  at  input  time 
Inductance  of  flux  tube  in  henrys  at  input  time 
Longitudinal  width  of  the  flux  tube  in  radians 

The  mean  drift  rate  for- electrons  in  the  flux  tube  at  input  time 
The  total  integrated  undecayed,  trapped  electrons  in  the  flux  tube 
at  input  time 

The  azimuthal  drift  current  in  amps  due  to  debris  at  input  time 
The  1st  time  derivative  of  ADI  in  amps/sec.  at  input  time. 


Output : 


T Output  time  in  seconds 

EL  L-shell  value  at  output  time 

Q Charge  on  the  flux  tube  in  coulombs  at  output  time 

QD  The  1st  time  derivative  of  Q at  output  time 

QDD  The  2nd  time  derivative  of  Q at  input  time 

F The  right  hand  side  of  equation  43 

DI  Drift  current  in  amps  due  to  both  debris  and  injected  electrons 

DID  The  1st  time  derivative  of  DI  at  input  time. 

Method: 

A simple  modified  Euler  numerical  method  is  used  to  solve  the  circuit  equation  for 
Q.  The  L-shell  value  is  updated  by  integrating  Equation  (42)  over  the  input  time 


Input : 

N The  number  of  values  in  the  input  array  D for  each  hemisphere 

ND  The  maximum  value  of  N and  the  1st  dimension  of  array  D 

D The  input  data  array  dimensioned  (ND,  2)  with  data  arranged  in 

each  hemisphere  from  cutoff  to  equator 

Output : 

R The  single  dimensioned  output  array  of  2N-1  values  where  R(l)  is  the 

southern  most  R(N)  is  the  equatorial  and  R(2N-1)  is  the  northern 
most  values. 

Method: 

N data  from  the  southern  hemisphere  of  D are  stored  sequentially  into  R and  N-l 

data  from  the  northern  hemisphere  of  D are  stored  in  reverse  order. 

t.  Function  PPLDC 


Purpose : 

This  routine  interpolates  to  find  the  mean  location  of  the  L-shell  value  of  the 
plasmapause  boundary  as  a function  of  local  time.  This  boundary  represents  a 
transition  between  open  and  closed  field  lines.  Within  this  boundary  the 
ionosphere  is  well  represented  by  a diffusive  equilibrium  model,  while  just 
outside  the  boundary  a convective  model  applies. 

Input : 

TL  The  local  time  in  hours  from  midnight 

Output : 

PPLOC  The  L-shell  boundary  of  the  plasmapause 

Method : 

The  plasmapause  position  is  interpolated  for  in  time  from  a table  obtained  from 
the  observations  of  D.  L.  Carpenter  (Reference  12). 
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u. 


Subroutine  SINMLT 


Purpose: 

This  subroutine  computes  the  absolute  value  of  the  sine  of  the  magnetic 
latitude  at  a point  given  by  the  ratio  of  the  magnetic  intensity  at  the  point 
to  the  magnetic  intensity  at  the  equator. 

Input: 

BTBEQ  The  ratio  B/Bo 

Output: 

SMLT  The  value  | sine  (magnetic  latitude)) 

Method: 

Subroutine  SINMLT  computes  the  absolute  value  of  x = sin  \m  given  the  ratio 
B/B0  of  magnetic  intensity  of  the  point  in  question  to  that  at  the  equator. 
Newton's  method  is  used  where  the  solution  to  the  equation 

f(x)  - 1 + 3x2  - (1-X2)6  (B/B0)2  - 0 (173) 

is  sought.  An  initial  guess  for  x=xj  is  taken  to  be  Xj  - [(  B/BD)2  + l] / 

[6(B/B0)2  + 3])1/2  if  B/BQ  < 2.66,  or  xl  - (1  - 3(B0/B) l>6) 1/2  if  B/Bq  ^2.66. 

Thus,  given  the  initial  value  of  x the  subroutine  iterates  using  the  relation- 
ship xi+1  = Xt  - f(xi)/f1(xi),  i-1,  2,  3 ...  until  |xi+1  - xjJ  1 10'6,  where 


v.  Subroutine  SLATSQ 

Purpose : 

This  subroutine  computes  the  square  of  the  dipole  latitude  corresponding  to 
the  mirror  point  of  an  electron  with  a given  equatorial  pitch  angle. 

Input : 

K • 0,  if  SLAT2  is  input 

• 1 , if  SLAT 2 is  not  input 
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SPA 

SPA2 

SLAT2 


The  sine  of  the  equatorial  pitch  angle 

The  square  of  SPA 

The  initial  guess  if  k = 0 


Output : 


SLAT  2 


Method: 


The  square  of  the  sine  of  the  dipole  latitude 


Subroutine  SLATSQ  computes  n 2 given  sin2a  by  an  iterative  method  using  the 


relationship 


sin2a0  = (l-nm2)3/^  + 3nm2)172. 


(174) 


This  method  consists  of  computing  a first  approximation  for  using  the 
equation, 

nm2i  = (23/16)-  V (23/16)2  - 4 (1231/1024)  (T/T0-l), 

2 (1231/1024) 


T <_  T*  = T 1 + 1/4  (1024)  (23)2 

(1231)  (16) 


1 , when  T > T* 


where 


T = 1.380173  - . 639693 (sin2a) 0,37 


(176) 


T0  - u/3/2  = 1.380173  - .639693. 


Equation  (174)  is  then  used  in  the  following  iterative  manner 


<i+l  - 1 - sin2a0(l  + 3^) 1/2  1/3 


to  obtain  a precise  value  of  r^2 
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w.  Function  SLFS 

Purpose: 

This  routine  computes  the  sine  of  the  dipole  latitude  corresponding  to  a given 
distance  along  a specified  field  line. 

Input : 

EL  The  L-shell  value  of  the  field  line 

S The  distance  in  km  to  the  equator  from  the  point  of  interest 

Output: 

SLFS  The  sine  (dipole  latitude) 

Method: 

An  initial  guess  is  provided  by  the  equation 

«0  * | S/LRe|  (179) 

This  guess  is  constrained  to  be  less  than  1.380173  and  then  refined  using  a 
quadratic  approximation 

a,  * aQ(l  - ,199579o0)  (180) 

This  is  followed  by  two  iterations  using  the  following  equation. 

i 

I 

«1+1  =*  j a1+(2a0-ln(*/3a1+R)//3)/R  (181) 

where  SLFS  * a and  R = /l+3«2  j 

x.  Subroutine  TANDHI 

Purpose : 

This  subroutine  computes  a mean  ionospheric  temperature  and  the  0+/H+  transition 
height  (the  altitude  where  n(0+)  ■*  n(H+))  as  a function  of  magnetic  latitude 
and  local  time. 
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Input: 

M 


DML 

TL 

Output: 

TIK 

HIKM 


= 1,  if  TIK  and  HIKM  are  both  to  be  computed 
= 2,  if  only  TIK  is  to  be  computed 
= 3,  if  only  HIKM  is  to  be  computed 
The  magnetic  latitude  in  deg. 

The  local  time  in  hours  from  midnight 


The  ion  temperature  in  deg.  K 
The  0+/H+  transition  altitude  in  km 


Externals: 


VOFY 


Method: 

The  data  used  in  this  code  is  taken  from  the  work  of  J.  E.  Titheridge 
(Reference  13).  The  ion  temperature  is  computed  for  daytime  (noon)  and 
nighttime  (midnight)  from  previously  stored  spline  interpolating  coefficients 
as  a function  of  the  absolute  value  of  magnetic  latitude.  A simple  linear 
interpolation  in  local  time  is  made  between  these  two  values  to  obtain  the 
temperature  at  the  input  local  time. 

The  0+/H+  transition  altitude  is  similarly  computed  as  a function  of  magnetic 
latitude,  however,  we  are  forced  to  assume  that  the  day  night  variation  which 
we  can  compute  for  the  northern  hemisphere  also  applies  in  the  southern 
hemisphere. 

y.  Test  Cases 

Three  test  cases  for  the  tube-motion  model.  Including  the  new  input  data  have 
been  provided.  These  test  cases  are  for  ARGUS-3,  STARFISH  and  TALL  BEAR.  All 
of  these  runs  use  AMBCAF  to  compute  the  tube  capacitance  and  resistance,  OEBCON 


i 


to  compute  the  tube  motion  and  debris  distribution  and  the  old  SPECTER 
INJECT  and  FLUXS  code  to  compute  drift  dilution  and  decay  of  the  injected 
electrons.  The  STARFISH  test  case  also  exercised  the  new  DSTRIB  code  to 
estimate  the  redistribution  of  the  injected  electrons  by  the  electric  field 
in  the  tube. 

The  input  to  AMBCAP  specifies  the  geographic  altitude,  latitude  and  longi- 
tude, the  local  time  and  the  day  of  the  year.  The  DEBCON  input  is  yield, 
kinetic-yield  fraction,  burst  radius,  NRL  reference  plane  distances, 
final  time  and  time  increment.  The  ARGUS-3  case  is  a low-yield,  high- 
latitude,  southern  hemisphere  run;  STARFISH  is  a medium-yield,  low- 
latitude,  northern  hemisphere  run,  and  TALL  BEAR  is  a high-yield,  high- 
latitude,  northern  hemisphere  test.  In  all  cases  the  tube  motion  time 
increment  was  chosen  so  that  the  details  of  the  debris  flux  and  velocity 
profiles  specified  in  the  reference  planes  could  be  resolved.  The  burst 
radii  were  taken  from  the  results  of  the  LMSC  code  that  computes  the 
reference  plane  parameters  and  cross-sectional  area. 

A test  case  has  also  been  provided  for  the  new  saturation  model  based 
of  flux  limiting  due  to  strong  diffusion.  In  this  test  the  injection 
rate  is  100  MT/sec  for  a duration  of  6 hours.  This  rate  establishes  a 
flux  which  exceeds  J*,  the  equilibrium  flux,  except  at  very  low  L values 
where  the  decay  rate  is  very  high  due  to  collisions.  The  flux  is  iso- 
tropic and  constant  along  the  magnetic  field  until  after  the  duration  of 
the  nuclear  exchange  (6  hours).  It  then  decays  at  a diminishing  rate, 
approaching  the  normal -mode  decay  rate  asymptotically,  as  discussed  in 
Section  IV. 4. 

The  results  of  the  test  cases  are  described  in  the  detailed  output  listings 
which  have  been  provided. 
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SECTION  III 


1 


COMPARISON  OF  COMPUTED  FLUXES  WITH  TESTS  DATA 

In  this  section  the  distributions  of  trapped  electrons  computed  for  the 
Argus-3  and  Starfish  bursts  are  described  and  compared  with  the  available 
data  on  the  radiation  belts  produced  by  those  nuclear  tests. 

1.  ARGUS -3 

During  the  Argus  tests,  three  identical  devices,  Argus  1,  2,  and  3,  were 
detonated  at  altitudes  greater  than  about  200  km  in  the  region  of  the  South 
Atlantic  geomagnetic  anomaly  in  late  summer,  1958.  The  yields  were  in  the 
range  1-2  KT.  Each  device  produced  an  intense,  narrow  radiation  belt  that 
persisted  for  several  weeks.  The  belts,  at  least  at  early  times,  were  quite 
similar. 

For  the  Agrus-3  test,  we  received  from  NRL  the  initial  distribution  of  the 
hot  ions  (the  ionized  debris  and  accompanying  energetic-air  ions)  at  the 
reference  planes  described  in  Section  II. 3a.  The  fission-fragment  flux 
and  the  total  ion  flux  at  the  reference  planes  are  shown  in  Figure  19. 

The  velocities  of  the  ions  in  the  reference  planes  were  initially  about 
510  km/ sec  but  rapidly  diminished  monotonically , qualitatively  in  the 
manner  of  the  fall-off  of  the  velocities  of  the  ions  for  the  Spartan  bursts 
shown  in  Figure  17,  reaching  a value  of  about  80  km/ sec  10  seconds  after  the 
burst.  The  effective  azimuthal  width  of  the  debris  tube  was  0.044115  rad, 
and  the  resistance  was  taken  to  be  less  than  the  critical  value  discussed 
in  II. 2b,  subsection  (b). 

The  motion  of  the  debris  tube,  computed  by  using  the  above  data  as  input  to 
the  new  injection  model,  is  shown  in  Figure  27.  There,  the  L value  of  the 
center  line  of  the  tube  is  plotted  against  time.  As  expected  for  small  dis- 
placements of  the  tube,  the  motion  is  oscillatory.  The  thickness  of  the 
intense  region  of  the  resulting  radiation  belt,  according  to  the  model,  can- 
not be  less  than  the  width  AL  of  the  debris  tube,  which  is  AL“L0AiJ>=  .10RE. 


2.30 


Figure  27.  Motion  of 


If  the  source  of  the  injected  electrons  remains  strong  in  the  tube  at  times 
greater  than  15  or  20  seconds,  the  displacement  of  the  tube  toward  higher 
L values  would  spread  the  electrons  over  a wider  range  in  L.  The  thickness 
in  L of  the  Argus-3  shell,  based  on  a measurement  by  the  Explorer  IV 
satellite  at  B=.23  gauss,  about  four  hours  after  the  burst,  was  about  0.146 
Rg  (Ref.  14).  This  thickness  is  consistent  with  the  width  of  the  debris 
tube  and  the  effective  duration  of  the  electron  source  during  the  tube 
motion. 

Unfortunately,  the  L grid  selected  to  read  out  the  fluxes  of  the  injected 
electrons  was  not  fine  enough  to  resolve  the  shell  thickness.  The  maximum 
flux  was  at  an  L value  near  burst  L value.  At  the  next  grid  joint,  dis- 
placed from  the  first  by  AL=.07Rg,  the  flux  had  dropped  off  substantially. 

The  flux  computation  in  the  belt  has  been  compared  only  with  a measurement 
made  with  Van  Allen's  Geiger  counter  (Channel  3)  during  the  Explorer  IV 
traversal  of  the  belt  4 hours  after  the  burst.  This  detector  had  a nominal 
threshold  energy  of  3 MeV  and  a geometric  factor  of  about  0.6  cm  . Its 
counting  rate  at  B-,23  gauss,  as  given  in  a graph  in  Ref.  14  depicting  the 

decay  rate  of  the  belt,  was  about  4x10^ /sec.  Hence,  the  measurement  implies 

4 -2  -1 

an  omnidirectional  flux  of  6x10  cm  *sec  of  electrons  of  energy  greater 

than  about  3 MeV.  The  computed  flux  in  the  most  intense  region  of  the 

belt  just  after  the  flux  became  uniform  in  longitude,  at  a B value  near 

7 2 

.23  gauss,  was  4.7x10  1 MeV-equivalent  electrons  per  cm  per  sec.  This 

flux  corresponds  to  a fission  beta  spectrum  of  electrons  of  energy  greater 
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than  3 MeV  of  about  4x10  cm  'sec  , an  order  of  magnitude  greater  than 
the  measured  value  given  above.  This  discrepancy  is  not  regarded  to  be 
excessive  in  view  of  the  uncertainties  in  the  measurement  and  in  the  fission- 
fragment  distribution.  Moreover,  the  comparison  was  made  at  a low  altitude 
where  the  flux  was  changing  rapidly  in  time  and  in  B value.  Further  com- 
parisons of  the  computed  fluxes  with  experimental  data  are  needed.  Such 
comparisons  will  be  facilitated  by  the  results  of  a separate  effort  funded 
by  the  AFWL  to  assemble  and  assess  all  of  the  data  obtained  in  the  high 
altitude  nuclear  tests. 


In  the  Argus  tests  the  debris  tube  moved  so  slowly  that  the  redistri- 
bution of  the  electrons  was  negligible.  Indeed,  all  of  the  available 
measurements  indicated  that  the  energy  spectra  of  the  electrons  were 
indistinguishable  from  the  equilibrium  fission  beta-decay  spectrum. 

2 . STARFISH 

The  fission-fragment  flux  in  the  reference  planes  for  the  Starfish  burst, 
determii.ed  from  the  data  provided  by  NRL,  is  shown  in  Figure  28.  The 
dashed  curve  in  that  figure  shows  the  ratio  of  the  fission-fragment  flux  to 
the  total,  energetic-ion  flux.  At  the  earliest  time,  about  0.2  sec,  at 
which  the  hot  ion  flux  appears  at  the  reference  plane,  the  flux  ratio  is 
.00917,  and  the  hot  ion  flux  consists  almost  entirely  of  debris  at  that 
time.  The  hot  air  predominates  at  later  times  as  mentioned  previously. 

The  velocity  of  the  hot  ions  at  the  reference  plane  as  a function  of  time 
is  shown  in  Figure  29.  The  azimuthal  width  of  the  debris  tube  was  .07778 
rad,  and  the  resistance  appropriate  for  the  undisturbed  ionosphere  along 
the  tube  motion  was  .774  ohms. 

The  debris  tube  motion  resulting  from  these  input  data  is  shown  in 
Figure  30.  Note  that  the  tube  accelerates  monotonically  as  it  convects 
through  the  magnetosphere  to  L values  beyond  the  trapping  limit.  This 
behavior  of  the  tube  motion  was  found  to  be  the  same  for  all  the  high- 
yield  (>1  MT)  bursts  for  which  input  data  was  received  from  NRL,  viz, 
for  Tall  Bear  and  the  Standard  Spartan  bursts  at  L-3  and  at  altitudes 
200,  300,  400,  and  600  km.  The  energy-differential,  omnidirectional 
fluxes  of  the  redistributed  electrons  at  various  L values  are  shown  in 
Figure  21.  As  discussed  previously,  the  corresponding  fluxes  that  are 
uniform  in  longitude,  neglecting  the  small  effects  due  to  decay  and  the 
"atmospheric-wiper"  action,  can  be  obtained  from  the  values  in  Figure  21 
by  multiplying  by  the  ratio  A«J> / 2it . Measurements  of  the  trapped-electron 
spectra  in  the  L range  1.25  to  1.70  were  made  by  West  et  al.  (Refs.  15  and 
16).  several  months  after  the  Starfish  test,  with  a mass  spectrometer  on 
the  Starad  satellite.  Their  spectra  at  various  L values,  together  with  the 
computed  spectra  (the  fluxes  of  Figure  21  multiplied  by  A4> / 2tt ) for  the  speci- 
fied L values,  are  shown  in  Figure  31.  Because  of  the  presence  of  natural 
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Figure  29.  Velocitv  of  Starfish  hot  ions  along  magnetic  field, 
at  reference  planes,  as  function  of  time. 
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Figure  31.  Experimental  spectra  of  radiation  belt 

electrons  following  Starfish  (Reference  15) 
The  solid-line  curves  without  the  open 
circles  are  the  computed  spectra  of  the 
redistributed  electrons  (see  Figure  21). 


trapped  electrons  and  the  decay  of  the  Starfish  electrons,  detailed  agree- 
ment is  not  expected.  However,  the  general  softening  of  the  computed 
spectra  at  increasing  L values  is  verified  by  these  data. 


1 


/ 


The  comparison  of  the  equatorial,  omnidirectional  flux  with  the  Telstar 
data,  however,  is  not  so  favorable.  Here,  the  electron  distributions 
computed  with  the  new  injection  model  and  with  the  old  jetting  model  are 
both  in  disagreement  with  the  Telstar  data  at  L values  above  the  peaks 
of  the  distributions.  In  making  the  comparison,  the  effect  of  the  different 
electron  spectra  given  by  the  two  models  was  removed  by  computing  the 
response  of  the  Telstar  detector  to  the  fluxes  by  integrating  the  counting 
efficiency  of  the  detector  given  in  Ref.  17  (channel  3,  with  lower  and 
upper  pulse-height  edges  at  390  and  615  keV,  respectively)  over  the  energy 
spectra  of  the  fluxes.  The  resulting  counting  rates  are  shown  in  Figure  32. 
The  open  circles  are  the  counting  rates  of  the  Telstar  detector,  at  the 
equator,  two  days  after  Starfish.  The  solid  curve  gives  the  counting 
rate  expected  for  the  flux  computed  with  the  new  injection  model.  Here, 
decay  is  not  included.  The  detector  response  was  computed  by  using  the 
spectra  shown  in  Figure  21  with  the  spectral  intensities  reduced  by  the 
factor  .07778/2w.  The  broken-line  curve  in  Figure  32  is  the  counting 
rate  expected  if  27%  of  the  debris  jets.  This  is  the  magnitude  of  the 
jetting  expected  at  early  times  when  the  ram  pressure  of  the  debris 
exceeds  the  magnetic-field  pressure.  The  broken-line  curve  was  computed 
using  the  old  jetting  model,  assuming  the  equilibrium  fission  beta-decay 
spectrum,  and  applying  a two-day  decay  factor.  The  data  should  be  compared 
with  the  sum  of  the  computed  counting  rates,  which  are  approximately  the 
rates  expected  if  the  debris  both  jets  at  early  times  and  is  carried  upward 
by  the  tube  motion  at  later  times.  As  shown  in  the  figure,  the  computed 
counting  rates  are  in  fairly  good  agreement  with  the  data  at  L<1.6.  How- 
ever, at  the  higher  L values,  the  Telstar  data  lie  higher  than  the  computed 
values  by  an  amount  that  increases  with  increasing  L value. 

It  is  well  known  (see  e.g.  Ref.  18  and  14)  that  the  fluxes  measured  by 
Telstar  were  much  higher  than  those  inferred  from  the  Injun  satellite  at 
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Figure  32.  Counting  rate  of  Telstar  detector  (open  circles)  near  equator 
as  function  of  L,  together  with  counting  rate  of  detector 
due  to  fluxes  computed  with  new  injection  model  (solid  line) 
and  jetting  model  (broken  line). 
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L>1.4.  It  was  thought  that  this  discrepancy  might  be  accounted  for 
by  the  redistribution  process  in  the  debris  tube  which  transports 
lower-energy  electrons  to  higher  L values.  But,  as  shown  in  Figure  32, 
this  redistribution  of  the  electrons  is  not  sufficient  to  account 
entirely  for  the  Telstar  measurements.  The  Telstar  data  have  not  been 
corrected  for  background  measurements;  the  natural  radiation  prior  to 
Starfish,  in  the  energy  sensitivity  region  of  the  detectors,  was  not 
known  since  Telstar  was  launched  a day  after  the  Starfish  burst. 

Vette's  AE-4  (1964)  steady-state  environment  for  E>500  keV  for  solar 
minimum  (Ref.  19)  indicates  that  less  than  one-tenth  of  the  detector 
response  can  be  attributed  to  the  steady  state  background  radiation. 

Of  course,  the  natural  radiation  in  the  region  of  the  slot,  where  the 
discrepancy  occurs,  is  highly  dependent  on  storm-time  activity;  the  flux 
can  easily  increase  by  an  order  of  magnitude  during  magnetic  storms. 
However,  no  major  magnetic  storms  occurred  for  several  months  prior  to 
Starfish.  The  discrepancy  might  be  explained  by  one  or  more  of  the 
following  processes: 

i)  the  injected  electrons  might  have  drifted  to  the  appropriately 
higher  L values  while  they  were  bunched  in  longitude,  as  dis- 
cussed in  Section  II. 4b,  or  even  at  later  times  due  to  the  fluting 
instability, 

ii)  the  oscillations  of  the  magnetosphere  due  to  the  explosion  might 
have  caused  the  trapped  electrons  in  the  outer  radiation  zone  to 
diffuse  into  the  slot  region,  and 

iii)  the  counting  rate  of  the  Telstar  detector  might  have  been  enhanced 
due  to  the  acceleration  of  ambient  electrons  by  processes  associated 
with  the  c . xosion. 

The  latter  possibility  was  first  suggested  by  Colgate  (Ref.  20);  he  esti- 
mated that  the  counting  rate  of  the  Telstar  detector  could  be  due  to  the 
acceleration  of  ambient  electrons  by  the  upward-moving  collision-free 
shock  produced  by  the  explosion.  Fapodopoulos , during  the  Specter  Review 
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held  at  NRL,  also  discussed  a shock  mechanism  capable  of  accelerating 
ambient  electrons  to  very  high  energies.  If  high  fluxes  of  such  elec- 
trons are  accelerated  to  energies  higher  than  about  3/4  MeV,  they  would 
be  of  Importance  to  the  SPECTER  program  and  should  be  assessed  carefully. 

The  combination  of  debris  jetting  at  early  times  (<.75  sec)  and  debris- 

tube  convection.  Including  redistribution,  appears  to  give  results  for 

Starfish  that  are  consistent  with  the  majority  of  the  experimental  data. 

The  modeling  of  these  processes  contains  the  correct  first-order  physics. 

The  jetting  is  expected  theoretically  and  it  has  been  observed  optically 

(Ref.  2fy;  it  also  explains  the  high  fluxes  (about  10^  effective  1 MeV 
2 

electrons  per  cm  per  sec)  at  low  altitudes  on  L=7  inferred  from  an 
analysis  of  the  response  of  the  heavily-shielded  Geiger  counter  on  the 
ARIAL  satellite  (Ref.  22)  using  the  calibrations  of  similar  detectors  by 
O'Brien  et  al.  (Ref.  23).  At  the  lover  L values,  the  results  are  in 
general  agreement  not  only  with  the  Telstar  data,  as  discussed  above, 
but  also  with  Injun  data.  Of  course,  further  comparisons  of  the  com- 
puted fluxes  with  the  data  are  needed  in  regions  off  the  equator. 

For  high-yield  bursts,  owing  to  the  equipotential-f ield-line  approximation 
used  in  the  debris-motion  model,  the  effective  inductance  of  the  field 
lines  and  the  expected  anomalous  resistivity  along  the  field  cause  the  poten- 
tial to  vary  along  the  magnetic  field.  The  effect  is  such  that,  at  early 
times,  while  the  majority  of  the  electrons  are  injected,  the  electric 
field  across  the  tube  at  the  lower  altitudes  is  not  as  high  as  the  model 
predicts.  Hence,  the  debris  density  and  the  attendant  electron  injection 
at  low  altitudes  are  somewhat  lower  at  high  latitudes  than  predicted  by  the 
model. 

The  hot-air  plasma  has  rendered  the  debris-tube  motion  for  high-yield  bursts, 
such  as  Starfish,  rather  insensitive  to  wide  variations  of  both  the  ionos- 
pheric resistance  and  the  tube  width.  Hence,  the  electron  distribution  given 
by  this  model  is  now  almost  entirely  dependent  on  the  plasma  data  provided 
by  NRL. 
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SECTION  IV 


THE  APPROACH  TO  SATURATION 

1.  THE  MECHANISMS  OF  SATURATION  IN  THE  RADIATION  BELTS 

Trapped  electrons  are  subject  to  the  influence  of  electromagnetic 
waves,  and  can  actually  cause  the  generation  of  waves  in  many  modes. 
Whenever  the  unstable  generation  of  waves  results  in  a loss  of  particle 
energy,  the  particles  are  redistributed  or  lost.  The  problem  of 
particle  loss  then  boils  down  to  a question  of  which  modes  of  wave- 
particle  interactions  dominate.  In  the  natural  radiation  belts  it 
seems  fairly  certain  that  a major  part  of  the  energetic  electron  loss 
is  due  to  interactions  with  waves  in  the  electromagnetic  cyclotron, 
or  whistler,  mode.  Several  mechanisms  have  been  proposed  for  the  way 
electrons  are  removed  by  interactions  with  waves.  In  the  Kennel  and 
Petschek.  theory  (Ref.  24)  self  interactions  limit  the  electron  fluxes: 
when  the  fluxes  increase  momentarily,  the  wave  generation  and  subsequent 
pitch  angle  diffusion  are  enhanced.  The  Lyons,  Thorne  and  Kennel  theory 
(Ref.  25  ) invokes  the  same  wave  modes,  but  the  waves  are  not  self- 
excited.  The  waves  are  supposed  to  fill  the  entire  trapping  region, 
and  can  originate  anywhere  within  the  magnetosphere.  More  recently, 
it  has  been  suggested  (Ref.  26  through  30  ) that  man-made  VLF  waves  can 
cause  a large  part  of  the  electron  losses  in  the  natural  radiation  belts. 
The  mechanism  is  superficially  similar  to  the  Lyons,  Thorne,  and  Kennel 
mechanism;  but  the  waves  may  be  coherent,  which  results  in  diffusion 
rates  significantly  higher  than  the  quasi-linear  rates  (Ref  .31  ). 

Though  there  is  some  controversy  about  the  true  importance  of  self- 
excitation  of  waves  in  the  natural  radiation  belts,  there  is  no  cuestion 
that  self-excited  waves  must  eventually  limit  the  total  number  of  electrons 
that  can  be  trapped.  It  also  appears  that  the  trapping  limit  must  occur 
long  before  the  particle  pressure  overwhelms  the  trapping  field  pressure 
(betat>>l).  In  fact  we  know  that  it  takes  a very  well  designed  experi- 
ment to  exceed  the  beta  = 1 limit  in  laboratory  plasmas.  The  aftermath 
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of  a high  altitude  nuclear  explosion  is  highly  chaotic,  and  very 
unlikely  to  provide  the  uniform  conditions  needed  for  very  high  beta 
ratios.  So,  it  is  worthwhile  to  ask  whether  the  Kennel  and  Petschek 
mechanism  limits  the  electron  fluxes  before  the  beta  criterion  is 
reached.  A secondary  concern,  to  which  we  have  no  answers  yet,  is 
whether  other  wave  modes  might  be  as  important  as  the  cyclotron  mode 
at  high  flux  levels.  (Electrostatic  modes  are  conceivable  candidates 
for  a flux  limiting  mechanism.) 

One  criterion  relevant  to  the  importance  of  various  wave  particle  inter- 
actions is  their  characteristic  electron  removal  times.  High  beta 
instabilities  generally  result  in  hydromagnetic  waves,  which  cannot 
remove  electrons  faster  than  an  Alfv4n  wave  can  travel  along  a field 
line.  Alfv£n  wave  travel  times  are  generally  of  order  several  minutes 
in  the  cuter  radiation  belts.  These  times  seemed  short  enough  that 
a beta  criterion  was  used  in  the  early  versions  of  Specter  to 
limit  the  electron  flux  at  saturation.  The  characteristic  time  for  cyclotron 
wave  interactions,  on  the  other  hand,  is  the  strong  diffusion  lifetime. 

At  saturation  there  is  a geometrically  determined  limit  to  the  time  taken 
for  a trapped  particle  to  diffuse  from  a pitch  angle  near  90°  to  a pitch 
angle  in  the  loss  cone.  In  strong  diffusion  the  lifetime  is  independent  of 
the  wave  amplitude.  The  electrons  cannot  be  removed  faster  than  the  strong 
diffusion  lifetime,  which  in  the  radiation  belts  is  about  10  to  100  seconds 
(Ref.  4).  This  lifetime  is  generally  somewhat  smaller  than  the  Alfv£n 
travel  time,  so  it  is  reasonable  to  suppose  that  the  Kennel  and  Petschek 
mechanism  does  indeed  limit  the  saturation  fluxes. 

A major  criticism  of  the  Kennel  and  Petschek  theory  is  that  it  requires 
containment  of  the  waves.  Whistler-mode  waves,  however,  often  propagate 
at  large  angles  to  the  magnetic  field.  It  is  not  clear,  therefore,  whether 
a thin  shell  of  artificially  trapped  electrons  would  be  subject  to  the  Ken- 
nel and  Petschek  limit.  The  application  of  a Kennel  and  Petschek  limit 
to  the  debris  tube,  before  drift  dilution  of  the  electrons,  is  also  some- 
what questionable.  It  might  be  supposed  that  the  debris  tube  provides  an 
effective  duct  for  whistlers,  but  the  wave  reflection  coefficient  and  many 
other  critical  parameters  remain  unknown. 
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The  aim  here  is  to  develop  a loss  model  wherein  the  saturation  electron 
loss  rate  depends  only  on  the  flux  level.  At  low  flux  levels  the  loss 
rate  must  revert  to  the  natural  decay  rate.  The  model  may  be  faulty  at 
very  early  times  — due  to  the  possibility  of  hydromagnetic  instabilities 
and  to  our  imperfect  understanding  of  wave  propagation  in  the  debris  tube 
— but  a single  injection  event  is  very  unlikely  to  reach  saturation. 

Where  the  new  model  becomes  particularly  useful  is  in  the  multi-burst  case, 
where  many  injection  events  can  occur  within  time  scales  of  minutes.  The 
total  field  energy  in  the  trapping  field  is  equivalent  to  about  150  MT  (Ref.  32). 
A relatively  small  number  of  high-yield  explosions  might  exceed  the  beta=l 
limit.  Whether  saturation  ever  does  occur  will  be  the  subject  of  the 
following  subsections. 


2.  THE  NON-LINEAR  MODEL  OF  SATURATION 


Schulz  (Ref.  33)  has  proposed  that  the  particle  flux  must  decay  as 


dJ  XJ 
dt  — 1+At 

s 


+ S 

o 


(182) 
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where  X is  a diffusion  loss  rate  (or  reciprocal  diffusion  lifetime)  and 
Tg  is  the  strong  diffusion  lifetime.  Eq.  (182)  is  applicable  to  the  Omni- 
directional flux  at  the  equator,  J.  The  first  term  on  the  right  contains 
the  effects  of  self-excitation  of  waves,  and  is  limited  at  high  diffusion 

rates  to  -J/t  . The  remaining  term,  S , contains  the  effects  of  particle 
s o 

sources;  in  equilibrium 


So 
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1+A  T 
o s 


(183) 


Typical  values  for  the  equilibrium  flux,  Jq  , and  equilibrium  (parasitic) 
diffusion  rate  are 


in4  -2  “I 

10  cm  sec  < 

- ..+6  -2  -1 

J <10  cm  sec 

(184) 
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10  ^sec  * < X 

in-5  -1 
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(185) 

i 


j 

i 


I 


\ 


121 


Kennel  and  Petsche(Ref.  24)  proposed  that  the  diffusion  rate  increases  with 
flux,  and  that  A is  proportional  to  the  wave  amplitude  squared.  The  differen- 
tial equation  for  A must  contain  terms  related  to  the  wave  growth  rate,  the 
wave  absorption  rate,  and  an  external  source  of  waves  to  account  for  parasitic 
diffusion.  Schulz  proposes  an  equation  for  A that  reduces  to 

£ i - Q<1')o)  <186> 


Though  the  model  is  an  extreme  simplification,  Eq.  (5)  probably  represents 
the  physics  fairly  well  to  first  order  in  J and  A.  The  parameters  P and  Q 
are  related  to  the  Kennel  and  Petschek.  limiting  flux,  J*,  the  growth  rate 
at  the  limit,  y* , the  wave  group  velocity  Vg,  and  the  wave  reflection  coef- 
ficient, R,  thus 
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y*  % .1  sec 


-1 


T*  a 1 ftUr  —2  —1 

J*  ik  10  L cm  sec 


Q % -V 
InR  % 


lnR  a 1 
8 lre  % 1 

-3 


(187a) 

(187b) 

(187c) 

(188a) 

(188b) 


A computation  of  the  group  velocity  for  waves  resonating  at  the  equator  with 
2 MeV  electrons  gives 


Q % .2/L  sec 


(189) 


Q/P  % 10nL"5 


(190) 
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The  actual  value  of  Q/P  is  not  very  critical;  it  suffices  to  most  purposes 
to  say  that  Q/P  is  of  order  J*. 

Before  proceeding  to  detailed  solutions  of  the  above  differential  equations, 
one  should  examine  the  behavior  of  the  solutions  near  equilibrium. 
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Let 
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X = (J  -Jq)  / Jo 

(191) 

y = ( A-  A ) / A 
o o 

(192) 

The  differential  equations  reduce  to 


iJL  = . = _ PJo  _ g x = -c  + C -3. 

dx  y+yx+x  1 2 y+yx+x 


(Note  that  the  actual  cases  to  be  solved  have  >>  C^>>  1.) 


(193) 


The  loci  of  constant  slope  y " = A are  the  curves 
C -C  -A 

yA  • -*1  <194> 

In  Figures  33  - 39  some  solutions  are  plotted  for  various  values  of  and 
C 2«  along  with  the  loci  of  constant  slope  (the  loci  of  constant  slope  are 
labeled  with  the  angular  slope  in  degrees). 


If  a trajectory  y(x)  passes  through  the  origin  its  slope  must  be  equal  to  the 
limit  of  y/x,  or 


dy 


_A 

dx 
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C +A 

w 


^2-Cl  C2~C1_^  2 1/2 

— Y 1 ±K— y~)  - c^1'2 


(195) 

(196) 


Eq.  (196)  is  complex  for 

c2  < (Z^+l)2 


(197) 


which  leads  to  periodic  solutions  that  spiral  forever  about  the  origin  x-y-O. 
Figures  33,  34,  and  37  show  examples  of  periodic  solutions.  These  solutions 
are  not  likely  in  the  realistic  situations  of  interest  here.  When  C2  increases 
sufficiently  A becomes  real,  with  two  values 
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Figure  38.  Flux  vs.  diffusion  rate  trajectories. 
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The  larger  slope  is  only  a limit  case  (the  near  vertical  curve  in 

Figure  30  is  an  example);  the  majority  of  solutions  enter  the  origin  with 

a slope  near  There  are  also  interesting  classes  of  trajectories  to 

the  left  of  the  limit  curve  that  overshoot  the  origin  and  approach  from  the 

left  side  (J<Jo,X<X  ). 

o 

Note,*‘too,  that  there  is  an  inflection  point  wherever 


dyA  2 

^ = A = [C2-C1-A]/(C2-C1-A-x(C1+A)] 


(199) 


The  locus  of  inflection  points  is  plotted  in  Figure  36  as  a dashed  line. 

The  inflection  points  on  Figure  39,  however,  lie  so  near  the  A=0  curve  that 
most  trajectories  turn  toward  the  zero  slope  curve  and  follow  it  closely 
for  the  major  part  of  their  descent.  These  are  the  trajectories  of  primary 
interest  in  the  present  context. 

Eqs.  (182)  and  (186),  with  their  widely  dissimilar  characteristic  time  scales, 
constitute  a pair  of  "stiff"  differential  equations.  An  excursion  from 
equilibrium  quickly  causes  the  growth  of  waves  within  several  wave  growth 
periods;  the  wave-particle  system  then  subsides  back  toward  equilibrium  at 
a rate  determined  by  diffusion.  The  second  part  of  the  J,X  trajectory  is  the 
only  one  that  concerns  us  here.  But,  from  be  preceding  discussion  we  deduce 
that,  if  P is  very  small,  the  decay  follows  very  closely  the  zero-slope  curve: 


dlnA 

dlnJ 


0 = 


PJ  X-(1+Xxs) [Q(X-Xo)+PSo]J 

-x2+s  (1+Xt  )X 
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or 


JX  ^ (1+Xt  )$(a-a  )+S  ] 
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(200) 
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One  may  substitute  J*  for  Q/P,  though  the  appropriate  constant  may  not 
be  exactly  identical  with  the  Kennel  and  Petschek  limit.  Finally,  one  has, 
for  X as  a function  of  J 
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(202) 
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*<•»  S - 27"  <1Ti>+hr  <147*)2+  r(1-3?»1/2 
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For  small  J: 

X(J)  ~ x0j^Tj“>  J <<  J*  (203) 

Near  J*: 


X(J)  % t—  (l-§  )]1/2,  J % J*  (204) 

T J " * 

S 

For  large  J 

X(J)  * i-  (^-1)  J » J*  (205) 

s 

Consider  now  the  extreme  Cases  of  strong  and  weak  diffusion.  In  weak 
diffusion,  where  J <<  J*,  and  X <<  1/t  , Eq*  (203)  predicts  exponential 
decay  at  a rate  only  slightly  different  from  the  natural  diffusion  rate 
X . The  strong  diffusion  case  is  described  by  Eq.  (205),  which  gives 

J(t)  'v.  J*+(Jo-J*)exp(-t/Tg)  (206 

The  flux  thus  decays  at  a nearly  exponential  rate  until  it  reaches  the  satura- 
tion level,  whereupon  there  is  an  abrupt  transition  to  weak  diffusion  at 
a rate  that  quickly  approaches  Xq. 

Eq.  (202)  is  plotted  in  Figure  40  for  several  values  of  the  parameters.  The 
transition  to  strong  diffusion  as  J is  increased  is  very  striking.  For  values 
of  J less  than  J*  the  natural,  or  parasitic,  diffusion  rate  can  always  be 
used  with  high  accuracy.  For  fluxes  slightly  greater  than  J*  the  diffusion 
quickly  reaches  the  strong  diffusion  rate.  However,  at  J*  the  diffusion 
rate  assumes  an  intermediate  value.  This  has  important  consequences  that 
will  be  discussed  below.  It  also  seems  to  imply  that  Q/P,  which  is  usually 
slightly  less  than  J*,  is  an  important  parameter  itself. 

Using  Eq.  (202)  it  is  a simple  matter  to  integrate  Eq.(182)  to  obtain  the 
time  history  of  the  flux  (it  is  not  a simple  matter  to  directly  integrate  the 
pair  of  differential  equations,  otherwise  much  of  the  above  discussion  would 
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have  been  unnecessary).  Figures  41  and  ^2  show  several  solutions  for  J(t);  the 
two  figures  are  identical  but  for  their  different  time  scales.  The  flux  does 
decay  at  nearly  the  strong  diffusion  rate  until  it  reaches  J*;  it  thereupon 
decays  at  the  natural  rate.  So  the  use  of  Eqs.  (182)  and  (202),  applied 
to  the  electron  fluxes  in  Specter,  leads  naturally  to  a saturation  value 
which  is  difficult  to  exceed  for  long  periods. 

3.  ATTAINABLE  SATURATION  LEVELS 


The  saturation  flux,  J*,  is  not  an  absolute  limit  that  can  never  be  exceeded. 

If  new  electrons  are  injected  at  a high  enough  rate,  a flux  much  larger  than 
J*  can  be  sustained  as  long  as  the  injection  source.  This  is  in  marked  con- 
trast with  the  former  Specter  saturation  model  where  an  absolute  limit  was 
assumed.  The  values  used  for  J*  are  based  on  the  best  available  data,  which, 
in  some  cases,  are  good  only  to  within  an  order  of  magnitude. 

To  inject  electrons  up  to  a level  J £ J*  it  is  only  necessary  to  overcome 

the  natural  diffusion  rate,  which,  for  MeV  electrons,  is  of  the  order 

X^  % 10  sec  The  time  scale  of  artificial  injection  events  is  always  less 

than  X but  usually  somewhat  greater  than  the  strong  diffusion  decay  rate,  1/t  . 
o s 

For  Tg  Schulz  (Ref.  4)  finds  values 


t ^ .0'.  L4/  J % .04  L4  (207) 

s C 

The  flux  from  a series  of  isolated  events  would  be  ecpected  to  build  up 
slowly  until  a flux  near  J*  is  attained. 

> 

To  maintain  (or  generate)  a saturated  flux,  J ^ J*,  it  is  necessary  to  set  up 
an  intense  electron  source,  S';  Eq.  (182)  then  becomes 


dt  - 


XJ 

1+Xt 
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+ S +S 
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(208) 


If  the  build  up  rate  does  not  exceed  the  wave  growth  rate  (of  order  .1  to 
1 sec  1) , the  diffusion  rate  will  continuously  adjust  to  the  values  given 
by  Eq.  (202).  Even  if  the  wave  growth  rate  is  exceeded  for  a short  time, 
the  diffusion  rate  will  quickly  relax  to  those  values  when  the  injection 
rate  is  reduced,  or  when  the  sustainable  flux,  J',  is  reached.  Eqs.  (182) 
and  (202)  give  for  J': 
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J ' % S'[J*+S'T  ]/[  X J*+S  '] 
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The  source.  S',  is  related  to  the  electron  injection  rate, 
injection 
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a cos  a d( 
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dN 

dt 


For  isotropic 


(211) 


where  N represents  the  total  number  of  electrons  trapped  in  a tube  of  unit 
cross  sectional  area  at  the  equator,  and  aQ  is  the  equatorial  pitch  angle. 
The  injection  rate  is  usually  known  in  units  of  MT  fission  g yield  per 
second.  If  W is  the  injection  rate  in  MT  per  AL  per  second  (assumed 
uniform  in  longitude) , the  relation  to  S ' is 


S'  a,  51.4  W x 6xl026  tt—TT  % 1‘2  x 1()10  T (212) 

^ 2nR  L 

L 

Properly  W should  contain  the  trapping  efficiency,  which  is  usually  less 
than  10%.  The  sustainable  flux  is  then 


J ' 


(J*+4.7x108L2W) 
(W+8.3x10“11l2x  J* 


(213) 


The  time  taken  to  reach  J ' is  approximately  x'fy  eva^uate^  at  J '• 


Some  values  are  given  in  Table  2 for  a case  L % 4,  X % 10  8sec  \ 

8 —2  —1  ® 

J*  % 5x10  cm  sec  . It  is  evident  from  the  table  (also  from  Figure  40)  that 

J " reaches  a broad  plateau  near  J*  long  before  t'  approaches  t . This  is 

in  marked  contrast  with  the  customary  presumption  that  strong  diffusion  must 

prevail  when  J ' is  near  J*.  The  presumption  of  strong  differences  at  J* 

has  led  to  the  misleading  conclusion  that  the  saturation  limit  is  unlikely 

ever  to  be  reached  (Ref  34).  That  deduction  was  based  on  a diffusion  rate 

X 'v  1/t  , which  corresponds  in  the  example  of  Table  2 to  injection  rates 
> 8 -1-1 

W *v  1 MT  (AL)  sec  . But  injection  rates  4 to  5 orders  of  magnitude  lower 
can  give  fluxes  within  90%  of  the  saturation  value.  Near-saturation  can 
be  attained  with  realistic  injection  rates  as  low  as  10  ^ MT  sec  ^ % .04 
MT  hr  \ This  corresponds  to  an  actual  fission  yield  >.4  MT  hr 
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Table  2 

Sustained  Flux  Levels  Related  to  Injection  Rate  at  L = 4 


W (MT  AL  1sec_1) 

T.  / -2  -lv 

J'  (cm  sec  ) 

t (sec) 

W x T (MT  t 

10'9 

7.5  x 105 

io6 

.01 

i io-8 

7.4  x IO6 

9.9  x 105 

.099 

1 10'7 

6.5  x IO7 

8.7  x IO5 

.87 

! 

3.0  x IO8 

4.0  x 105 

4.0 

io-5 

j 

4.7  x IO8 

6.0  x ±0* 

6.0 

-4 

8 

10 

5.0  x 10 

3.2  x 10 

3.2 

io"3 

5.1  x IO8 

500 

5.0 

io-2 

Q 

5.8  x 10 

72 

7.2 

O 

1 

1.25  x IO9 

17 

17 

1 1 

4.0  x IO9 

11 

110 

j io 

7.6  x IO10 

10 

IO3 

2 

11 

4 

| 10 

7.5  x 10 

10 

10 

1 103 

12 

. 7.5  x 10 

10 

IO5 

J 


138 


4.  INCORPORATION  OF  SATURATION  EFFECTS  IN  SPECTER 


/ 


There  are  two  needs  for  a new  saturatioan  model  in  Specter.  One  is  in  the 
saturation  model,  where  saturation  is  assumed  to  prevail  for  a specified 
duration.  In  that  case  one  need  only  replace  the  former  beta  criterion  by 
the  new  criterion,  Eq.  (213).  The  adjustable  parameters  now  are  W and  the 
duration  of  the  saturation  event.  These  paraemters  have  readily  apparent 
physical  significance.  It  is  even  conceivable  that  W be  allowed  to  vary 
with  time,  in  which  event  J'  will  adjust  nearly  instantaneously  to  the 
changes . 


The  other  need  for  a saturation  model  is  in  the  decay  scheme,  where  the  decay 
rates  must  be  adjusted  continuously  when  the  flux  approaches  J*.  One  way 
this  can  be  achieved  is  to  carry  an  explicit  integration  of  Eqs.  (182)  and 
(202)  in  the  loss  model.  There  are,  however,  simple  approaches  that  give 
the  same  limiting  behavior  at  weak  and  strong  diffusion.  The  transition 
regime  is  not  well  enough  understood  that  the  Schulz  model  is  to  be  preferred 
over  all  others.  It  is,  though,  a valuable  guide  to  our  understanding. 

The  time  behavior  is  actually  represented  quite  well  by  the  formula 


J(t)  % J(0) [l-C*]exp(-t/x  )+J(0)C*exp(-t/t  ) 

S L 


j* 


C*  = [constant] 


J(0)+J* 


tt  = 1/A 

L O 


(214a) 

(214b) 

(215) 


That  this  form  of  decay  is  a good  representation  in  the  limits  is  readily 
apparent  from  Figures  41  and  42.  In  fact,  with  [constant]  % 1,  Eq.  (214) 

Table  3 lists  several  comparisons 
At  smaller  fluxes  the  loss 


is  a very  close  approximation  for  all  J > J*. 
for  the  case  J(0)  ^ 10^, 

time,  T=J/^jT*  approaches  more  quickly  than  in  the  explicit  solutions; 
but  either  representation  is  adequate  within  the  bounds  of  our  ignorance 


J*=108,t  =10,t=106. 

S la 


of  the  actual  physical  processes. 


The  decay  rates  of  Eq.  (214)  are  incorporated  in  the  revised  Specter  elec- 
tron decay  model.  The  coating  changes  are  minor,  and  affect  only  one  sub- 
routine — the  one  used  to  compute  less  rates. 
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Table  3 

The  Decay  of  a Saturated  Radiation  Belt 

t(sec) 


>?.c) 

J (cm 

-2 
j 

-1* 
sec  ) 

from  Eq. (182) 

from  Eq. 

-(214) 

0 

io10 

10.09 

10.10 

20 

1.43 

X 

109 

10.66 

10.7$ 

4o 

2.71 

X 

io8 

14.91 

15.86 

60 

1.14 

X 

108 

46.27 

83.11 

80 

9.24 

X 

107 

278 

7.59  x 

iou 

100 

8.96 

X 

IO7 

1.99  x IO3 

1.05  x 

IO5 

120 

8.92 

X 

IO7 

h 

1.44  x 10 

1.09  x 

105 

150 

8.92 

X 

IO7 

2.27  x IO5 

1.09  x 

105 

200 

8.91 

X 

IO7 

9.78  x IO5 

1.09  x 

IO5 
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SECTION  V 


INTERACTION  REGIONS  OF  MULTIBURSTS 

An  investigation  has  been  conducted  to  determine  under  what  circumstances 
the  simple  superposition  of  injected  fluxes  due  to  multiple  bursts  may  be 
applicable.  Specifically,  it  calls  for  the  boundaries  in  space  and  time 
that  define. 


i)  the  interaction  region,  wherein  the  detonations  are  not  sufficiently 
close  together  to  be  considered  as  a single  detonation  but  are 
sufficiently  close  together  to  affect  each  other's  injection  pro- 
cesses, 

ii)  the  summing  region,  wherein  the  detonations  are  sufficiently  close 
together  in  space  and  time  to  be  considered  as  a single  detonation; 


iii)  the  noninteracting  region,  wherein  the  electron  injections  are 
independent  and  unaffected  by  the  additional  detonations. 

First,  we  will  try  to  define  regions  (i)  and  (ii) . Region  (iii)  will  then 
be  specified  as  the  remaining  space  and  time  conditions  outside  of  (i)  and 
(iii). 


At  the  outset.  It  should  be  emphasised  that  a single  high-yield  burst  greatly 
disturbes the  atmrsphere-lonosphere-magnetosphere  system:  the  atmospheric 
heave  Inc rear  neutral  density  in  the  upper  atmosphere  by  several  orders 

of  magnitude  v hundreds  of  seconds;  the  heated  ionospheric  plasma 

fills  the  loca  tic  tube  of  force,  creating  an  anomously  high  plasma 

dansity  in  the  tube  that  persists  for  hours;  shock  waves  are  produced  that 
travel  several  times  around  the  earth;  the  entire  magnetosphere  oscillates, 
producing  a wida-band  of  micropulsationa  throughout  the  world;  and  even  the 
high-energy  protons  (55  MeV)  in  the  radiation  belt  are  redistributed  (see 
a.g.  Kef.  2).  In  the  stricteat  sense,  therefore,  the  interaction  region 
ahould  include  almost  all  the  epace-tine  acenarlos  envisaged  for  a nuclear 
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exchange.  In  this  study  we  will  attempt  to  include  only  first-order 
effects. 

1.  INTERACTION  REGION 

a.  Expansion  of  Collision-Dominated  Region. 

For  an  explosion  below  the  exosphere,  in  the  collision-dominated  region,  the 
expansion  of  the  nuclear  debris  will  be  stopped  within  a sphere  that  contains 
a mass  of  air  approximately  equal  to  the  mass  of  the  debris.  Owing  to  the 
buoyancy  of  the  fireball,  the  debris  will  also  be  accelerated  upward;  if  the 
fireball  radius  is  at  least  comparable  with  the  scale  height  of  the  atmosphere, 
the  fireball  will  undergo  a ballistic  rise  which  may  carry  the  debris  beyond 
the  exosphere.  However,  for  detonations  away  from  the  South  Atlantic  anomaly, 
very  few  fission  fragments  are  expected  to  reach  altitudes  above  about  1500  km, 
where  the  electron  trapping  efficiency  becomes  finite.  Therefore,  the  flux 
of  trapped  electrons  resulting  from  the  direct  injection  process  should  be 
negligible.  Nevertheless,  the  Teak  and  Orange  tests  did  establish  radiation 
belts.  It  is  believed  that  these  belts  resulted  from  pitch-angle  diffusion 
of  the  electrons,  due  to  a plasma  instability,  after  the  electrons  emitted 
by  the  fission  fragments  penetrated  the  residual  atmosphere.  Since,  at  the 
present  time,  the  Specter  codes  compute  trapped-electron  fluxes  due  to  the 
direct  injection  process  alone,  explosions  below  an  altitude  of  200  km  are  not 
included  in  this  study. 

In  a multiburst  situation,  however,  the  altitude  of  the  exosphere  is  dependent 
on  the  timing  of  the  nuclear  bursts  as  well  as  the  magnitude  of  the  attendant 
energy  deposition  in  the  atmosphere.  Data  on  the  increase  in  the  density  of 
the  atmosphere  at  high  altitudes,  in  time  and  space,  due  to  Standard  Spartan 
bursts  at  altitudes  of  175,  200,  225,  and  250  km  have  been  generated  by  NRL 
MRHYDE  code.  The  increase  in  the  altitude  of  the  region  where  collisions 
become  important,  as  a function  of  the  time  after  a 200-km  Spartan  burst, 
is  shown  in  Figure  43.  Note  that  the  altitude  of  the  collision-dominated 
region,  under  the  burst,  increases  substantially,  and  it  remains  high  for 
several  hundred  seconds.  The  collision-dominated  region  is  assumed  to  be  at  alti- 
tudes below  the  200-km  density  contour.  The  enhanced-density  region  is 
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Figure  43.  Altitude  of  200-km  atmospheric-density  contour  versus 
time  after  Standard-Spartan  burst. 
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bell-shaped  about  a vertical  axis  through  the  burst.  Since  the  motion  of 
the  debris  resulting  from  an  explosion  near  the  exosphere  is  not  known 
sufficiently  well  to  warrant  a detailed  specification  of  the  collision- 
dominated  region,  the  region  will  be  described  by  simplified  "square" 
functions.  Hence,  the  increased  height  of  the  region  will  be  taken  to  be 
Ah,  half  the  maximum  altitude  reached  by  the  200-km  density  contour  shown 
in  Figure  43;  the  diameter  of  the  region  will  be  taken  to  be  2R,  the  dia- 
meter of  the  200-km  density  surface  at  the  altitude  Ah;  and  the  time  inter- 
val of  the  enhanced  density  will  be  taken  to  be  At,  the  time  during  which 
the  altitude  of  the  contour  shown  in  Figure  43  exceeds  Ah. 


The  data  for  the  Standard  Spartan  bursts  at  different  altitudes  reveal  that 
the  intervals  At  are  quite  similar,  beginning  at  about  60  sec  and  ending  at 
about  480  sec.  However,  the  heights  Ah  decrease,  and  the  diameters  2R 
increase,  as  the  burst  altitude  increases.  These  variations  are  shown  in 
Figure  44 . 

For  bursts  of  different  yields,  but  at  a given  altitude,  the  following  reason- 
able assumptions  will  be  made: 

2 

i)  The  enhanced-density  volumes,  R Ah,  are  proportional  to  the 
yield.  Hence, 


R2Ah  Y 


(216) 


\ 


where  Y is  the  Standard  Spartan  yield  and  R and  Ah 
s s s 

are  the  dimensions  of  the  enhanced-density  region  given  in 
Figure  44 . 


ii)  The  shape  of  the  perturbation,  Ah/R,  remains  constant*  i.e.,  the 

atmospheric  expansion  pattern  for  a fixed  illumination  geometry 
remains  the  same.  Hence, 


Ah 

R 


(217) 
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Figure  44.  Effective  height-increase,  Ah,  and  diameter,  2R, 

of  collision-dominated  region  vs.  height  of  Standard' 
Spartan  burst. 
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The  above  equations  give, 
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i 


i 


i 


(218) 


(219) 

For  a burst  at  the  altitude  H*  of  yield  equal  to  5 MT,  for  example, 

Ah  ■ 110  km  and  R * 195  km.  If  a second  burst  is  detonated  in  this 
region  at  times  after  the  first  burst  in  the  range  60<At<480  sec, 
the  injection  efficiency  of  the  second  burst  will  be  impaired. 

If  the  disturbance  produced  by  a second  burst  overlaps  that  produced 
by  the  first,  an  enhancement  of  the  altitude  of  the  overlapping  region 
should  occur.  This  situation  is  depicted  in  Figure  45.  The  heights  and 

•ff  ^ 

yields  of  the  two  bursts  are  , Y^  and  Hj  , Y2;  the  parameters  of 

the  enhanced-denslty  region  produced  by  each  burst  alone  are  denoted  by 
the  subscripts,  1 and  2.  Ah^  is  the  height  of  the  enhanced-density 
region  over  the  overlapping  areas. 

The  altitude  Ah^  can  be  estimated  as  follows.  From  i) 

Rn2Ahn  " k(Hn)Yn>  n = X»2  (220) 


1/3 


Ah  - Aha(i-) 


and 


1/3 


R - Rs(f-> 


where,  as  indicated,  the  proportionality  constant  k is  a function  of 

* 

the  burst  altitude  H^.  The  average  value  of  k for  the  two  detona- 
tions is,  therefore. 


ave 


1 .RlAhl 

2 1 Y, 


R?Ah- 
+ -T-4) 


(221) 


atmospheric  density  produced  by  two  bursts. 


Thus,  for  the  combined  bursts, 


R,Ah,  = k (Y1  + Y0) 
3 3 ave  1 2 


(222  ) 


Similarly,  from  ii). 


Ah  Ah  Ah 

— m — ( + —) 

R3  2 V Rl  T R2; 


(223) 


Using  these  equations,  a solution  for  Ah^  can  be  expressed  as, 


Ah  ■ i (f  + f )^^  (Y  + Y )^^ 
s 2 '•12  21;  ' 1 2} 


where. 


Ah  R 

— (Ahx  + ^ Ah2) 


(225  ) 


and  f2^  is  obtained  from  (225)  by  interchanging  the  subscripts  1 and  2. 

* 

As  an  example,  for  Y^  = 5 MT  and  = 218  km,  Ah^  * 110  km  and  R^  * 195  km; 
for  Y2  = 10  MT  and  = 240  km,  Ah2  = 115  km  and  R2  » 251  km.  Hence,  if 
these  detonations  are  separated  in  time  by  less  than  480  sec,  then  in  the 
overlapping  spatial  region,  Ah^  * 145  km. 

b.  Region  of  Longitudinally-Bunched  Electrons 

An  impairment  of  the  injection  efficiency  also  occurs  if  a weapon  is 
exploded  in  the  region  of  a drifting  beta  tube;  i.e.,  in  that  region  where 
the  electrons  Injected  by  a previous  burst  are  still  highly  bunched  in 
longitude.  Such  a region  contains  an  eastward  electric  field,  which  tends 
to  redistribute  the  electrons  to  higher  L values,  and  magnetic-field- 
aligned  currents,  which  tend  to  neutralize  excess  charges.  A burst  in  that 


region  modifies  both  the  eastward-electric  field  and  the  current  patterns, 
thereby  altering  the  redistribution  of  the  electrons.  Furthermore,  the 
existing  electric  field  and  currents  modify  the  cross-L  motion  of  the 
debris  tube.  The  electron  loss  rate  due  to  wave-particle  interactions 
is  also  enhanced  because  of  the  increased  electron  flux. 

Initially,  the  azimuthal  sector  of  the  beta  tube  is  equal  to  the  width 

of  the  debris  tube,  A<j>,  used  in  the  injection  model.  At  later  times, 

owing  to  the  energy  spread  of  the  electrons,  the  width  of  the  sector 

increases  as  the  electrons  drift  toward  the  east.  We  will  regard  (i) 

the  center  of  the  bunched  electrons  to  be  at  the  magnetic  longitude,  6q, 

where  half  of  the  injected  electrons  are  at  d»4>  , (ii)  the  western  edge 

— o 

of  the  bunch  to  be  at  <j>^,  where  three-fourths  of  the  injected  electrons 
are  at  and  (iii)  the  eastern  edge  of  the  bunch  to  be  at  <J>2*  where 

one-fourth  of  the  injected  electrons  are  at  * Hence,  one-half  of 
the  electrons  will  be  within  the  magnetic  longitudes  and  ij^.  If 
N(L,t^)AL  is  the  total  number  of  electrons  injected  at  L in  the  interval 
AL  and  at  the  time  t^,  then  the  number  of  electrons  that  drift  past  <j>^  at 
the  time  t is 

oo 

N(^>_<(i1,L,t)AL  = N(L,t^)AL  f G(w)dw  (226  ) 


where  G(w)  is  the  energy  spectrum  of  the  redistributed  electrons  at  the 
eastern  edge  of  the  debris  tube.  It  is  given  by  the  sum 

exp(-w/w  ) 

G(w)  - (227) 

J j 

The  lower  limit  in  (226)  is  the  energy  of  an  electron  that  drifts  eastward 
to  i at  the  time  t;  i.e.,  it  satisfies  the  equation 


lD(Wl,L)(t-tL) 


(228) 
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where  <p  is  the  longitude  of  the  burst,  and  ^(v^,L)  Is  the  mean  azimuthal 
drift  velocity  at  L of  an  electron  of  energy  w^.  The  integration  in  (227) 
is  just  Skj  exP  (-Wj/Wj).  Hence.  (226)  can  be  solved  for  and 
corresponding  to  ratios  N($>$^,L,t)/N(L,tL)  of  3/4  and  1/4  respectively. 


The  interaction  region  can  then  be  specified  as. 


>2  - AL,  t - tL 


(229) 


where 


►l  ■ 4>  + + lD(w1,L)(t-tL) 


(230) 


>2  * ♦ + + ^j)(w2  ,L)  (t”t^) 


(231) 


t-f  < ff/2 

L - <L(w„l 


(232) 


c.  Region  of  Drifting  Debris  Tube 

Another  interaction  situation  occurs  when  a second  weapon  is  exploded  in 
a drifting  debris  tube  produced  by  an  earlier  burst,  or  when  any  debris 
tube  is  overtaken  by  another.  Such  a situation  also  leads  to  a disturb- 
ance in  the  electric  field  and  magnet ic-f ield-aligned-current  systems, 
hence,  to  a modification  of  the  injected -electron  distribution.  Again, 
owing  to  the  enhanced  electron  flux,  the  loss  rate  also  increases. 


The  Interaction  region  due  to  a burst  in  a drifting  debris  tube  may  be 
specified  as 


$ , (A^  +A$2)/2;  L,  (ALx  +AL2)/2;  t-t  , At 


(233) 


ft  Ik  ft 

where  $ , L , t are  the  coordinates  and  time  of  the  first  burst. 
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A4>  and  AL  are  the  dimensions  of  the  debris  tubes;  the  subscripts  1 and  2 
denote  the  parameters  due  to  the  first  and  second  bursts  respectively 


* 

L = L 


* 1 
- t )dt 


AL  - — 51  A*,  n = 1,2 
n L u 
o 

n 


AL,  + AL~ 

and  At  = — 

L(t-t  ) 


(234) 


(235) 


(236) 


L is  the  outward  drift,  dL/dt,  of  the  first  debris  tube. 

The  interaction  region  due  to  one  debris  tube  overtaking  another  can  be 
described  by  the  equations: 


, (A<)>1  + A<j)2)/2;  Lj^  < L < L2,  At  £ T2  - Tx 


where 


f m 

r = / ^ , n - 1,2 

n J Ln  ’ 


and  10 


T1  - T2  ; L1  - L2 


(237) 


(238) 


(239) 


Here,  L is  the  maximum  L value  for  persistent  trapping.  Computer  runs 
m 

are  required  to  determine  the  values  of  T for  various  initial  values  of  L. 
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2.  THE  SUMMING  REGION 

It  appears  that  the  summing  region  is  very  limited.  The  yields  of  two 
or  more  bursts  may  be  summed  only  if  one  or  more  bursts  occur  within  the 
region  where  the  debris  due  to  a previous  bur  t is  expanding  radially 
against  the  magnetic  field.  Hence,  the  bursts  must  be  separated  less 
than  about  a second  in  time  and  less  than  about  a single-burst -magnetic- 
bubble  diameter  in  space.  Furthermore,  the  spatial  separation  and  timing 
of  the  bursts  must  be  such  as  to  avoid  fratricide.  Here,  the  principal 
concerns  may  be  the  X-ray  fluence,  which  may  collapse  the  structures  within 
about  a millisecond,  and  the  neutron  flux,  which  may  substantially  reduce 
the  yield  of  a second  burst  due  to  overinitiation.  Since  the  travel  time 
of  the  neutrons  is  also  of  the  order  of  a millisecond  for  a magnetic- 
bubble-diameter  separation,  the  bursts  must  occur  in  the  restricted  region 
in  space  and  in  less  than  a millisecond  in  time  to  satisfy  the  summary 
condition.  It  seems  very  unlikely  that  such  conditions  can  be  met. 

3.  THE  NON-INTERACTION  REGION 

The  noninteraction  region,  to  first  order,  is  that  region  (Ax, At)  out- 
side the  regions  described  above  in  Sections  V.l  and  V.2. 
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Section  VI 


RADIAL  DIFFUSION  OF  ARTIFICIALLY  TRAPPED  ELECTRONS 
1.  THE  NEED  FOR  A NEW  TREATMENT  OF  RADIAL  DIFFUSION 

The  electron  loss  model  in  Specter  was  originally  based  on  the  theory  and 
observations  of  pitch  angle  diffusion.  Pitch  angle  diffusion  is  expressed 
in  a loss  rate,  which  represents  the  rate  at  which  trapped  electrons  are 
precipitated  into  the  atmosphere.  The  other  forms  of  wave-particle  inter- 
actions, that  lead  to  redistribution  rather  than  immediate  loss,  were  ig- 
nored. The  principal  process  that  redistributes  electrons  between  L-shells 
is  radial,  or  cross-L,  diffusion.  In  the  natural  radiation  belts,  radial 
diffusion  is  often  important,  but  it  seldom  dominates  the  loss  processes. 

Until  now  the  effects  of  radial  diffusion  have  not  seemed  important  enough 
to  include  in  Specter,  mainly  because  of  the  uncertainties  in  the  injec- 
tion and  initial  debris  distribution.  However,  with  the  development  of  an 
improved  injection  model  a new  consideration  of  radial  diffusion  is  called 
for. 

Radial  diffusion  has  usually  been  treated  in  an  approximation  where  pitch- 
angle  diffusion  is  accounted  for  by  a simple  loss  rate  (Ref.  35).  This 
leads  to  a diffusion  equation  for  the  flux  as  a function  of  L (35,36).  But 
the  first  adiabatic  invariant  must  be  assumed  constant  in  such  a formulation 
if  the  variation  of  energy  with  L is  to  be  computed.  This  contradicts  the 
known  fact  that  the  first  adiabatic  invariant  is  rapidly  altered  by  pitch 
angle  diffusion.  It  therefore  appears  that  a thorough  treatment  of  diffu- 
sion must  include  both  pitch-angle  diffusion  and  radial  diffusion  on  an  equal 
basis. 

In  the  case  of  artificial  radiation  belts  there  are  some  special  reasons  why 
radial  and  pitch  angle  diffusion  must  be  treated  simultaneously.  The  initial 
injected  distribution  is  a highly  non-equilibrium  distribution,  subject  to 
Initial  loss  rates  that  may  be  much  higher  than  the  late  time  rates  (Ref.  2,3). 
The  simplifying  assumption  of  a constant  loss  rate  (Ref.  35)  is  not  valid. 

A more  subtle  defect  in  the  early  treatments  has  to  do  with  the  dependence 
of  the  diffusion  rates  on  pitch  angle.  Pitch-angle  diffusion  — driven  by 
doppler-shifted  resonances  with  cyclotron  waves  — tends  to  favor  particles 


with  small  pitch  angles  (Ref.  25).  Electrons  with  pitch  angles  greater  than 
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60°  to  80°  do  not  respond  rapidly  to  pitch-angle  diffusion.  Radial  dif- 
fusion, on  the  other  hand,  favors  electrons  with  pitch  angles  near  90® 

(Ref.  37).  So  there  is  a possibility  that  radial  diffusion  could  move  a "pan- 
cake" distribution  near  the  equator  inward  or  outward  with  rates  higher  than 
predicted  by  the  simplified  calculations  with  constant  loss  rates,  indepen- 
dent of  pitch  angle. 


The  generalized  diffusion  equations  for  three  degrees  of  freedom  is  (Ref.  6) 
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where  F(M,J,f)  is  a phase  space  distribution  function  in  the  three  adiabatic 
invariants,  M.  J,  and  4>;  and  ^is  a Jacobian  relating  the  adiabatic  invar- 
iants to  three  coordinates  x^.  Walt  (Ref.  39)  suggested  that  one  of  the 
degrees  of  freedom  in  (240)  could  be  eliminated  by  choosing  to  be  more 
approximately  conserved  in  both  radial  and  pitch-angle  diffusion.  A use- 
ful variable  is 
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(241) 


This  variable  is  conserved  to  a high  degree  during  pitch-angle  diffusion, 
and  nearly  conserved  during  radial  diffusion.  In  radial  diffusion  the 
exponent  is  exactly  3 for  90®  pitch  angles  and  approaches  2 for  small  pitch 
angles  (Ref.  40).  For  most  radiation  belt  particles,  an  appropriate  n is  of 
order  2.6  to  3. 


The  use  of  an  assumed  constant  n is  subject  to  the  criticism  that  the  energy 
variations  are  not  precisely  accounted  for  if  the  diffusion  equation  is 
set  up  with  n constant.  This  is  not  necessarily  a serious  objection  for 
trapped  electrons  which  diffuse  predominantly  toward  lower  pitch  angles 
with  little  energy  change.  Only  in  narrowly  defined  regions  are  pitch- 
angle  and  radial  diffusion  equally  important;  usually  one  or  the  other 
clearly  predominates.  Only  if  a particle  can  pa.=s  back  and  forth  many  times 
between  the  regions,  can  the  approximation  of  n constant  have  serious 
consequences.  Uusually,  however,  once  an  electron  leaves  the  radial- 
diffusion  region,  it  is  thereafter  primarily  subject  to  pitch-angle 


i 
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diffusion.  Figure  46  shows  a concentrated  electron  distribution  with 
the  regions  schematically  indicated  where  each  mode  of  diffusion  dominates. 

A diffusion  trajectory  is  sketched  to  illustrate  that  the  diffusion  proces- 
ses are  usually  quite  independent. 

To  construct  the  diffusion  equation  in  L and  y =cos  a,  it  is  first  necessary 
to  obtain  the  Jacobian,^-.  The  first  and  second  and  third  adiabatic 
invariants  as  functions  of  y,  I.,  and  n,  are,  respectively 


« ■ Si'  n L3'n  a-,2> 
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(242) 
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2l  (244) 

$ = 2irBgRg/L 

where  B„  is  the  earth's  surface  field  at  the  equator,  and  s is  the 
(dimensionless)  distance  along  a field  line.  The  Jacobian  is  now  (Ref.  37) 
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where  T represents  the  quarter  bounce  integral. 
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* 1.3802  - .6397 (1-y) 


The  diffusion  operators  are  found  by  applying  the  chain  rule  to  the  partial 
differentials : 


/3  \ / 3 \ , (3-n) (1-y  J /3  \ 
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(247) 


The  customary  notation  (— )^  means  a partial  derivative  with  y held 

constant.  At  this  point  it  is  clear  that  n»3  is  an  especially  felicitous 

32  g2 

choice,  since  the  cross  derivatives,  and  , are  thereby  eliminated. 

Even  for  n less  than  3 the  cross  derivatives  are  small.  Ignoring  the 

cross  derivatives  leads  to  the  diffusion  equation 
3f  _ 1 fir  _ 3 f 1 
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(248) 


The  distribution  function  here  is  the  phase  space  distribution  at  the 
equator.  In  terms  of  the  directional  flux  j 
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Figure  46.  A typical  artificial  electron  distribution 
(contours)  in  L and  cosine  of  pitch-angle. 

The  dotted  lines  separate  regions  where 
radial  diffusion,  pitch  angle  diffusion, 
and  atmospheric  loss,  respectively,  dominate. 
A hypothetical  diffusion  trajectory  is  shown 
as  a dashed  curve . 
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f - j(E,y,L)/p" 


(249) 


(The  diffusion  coefficients  are  written  with  a bar,  D,  to  signify 
that  these  are  bounce  averaged.) 


It  is  curious  that  for  n = 8/3  the  L diffusion  operator  reduces  to  the 
earlier  form  (Refs.  1,2).  This  corresponds  to  equatorial  pitch  angles 
near  40°.  The  diffusion  operator  for  n=3  is  (Ref.  39) 


5/2  3 


LJ/i  «_|d  l-5/2  — I (250) 

3L  l LL  3L  J 

n 

Various  values  have  been  published  for  the  diffusion  coefficients.  For 
D we  have  used  a value  adjusted  to  match  the  values  of  Lyons,  Thorne, 
Kennel  (Ref.  25): 

°yy  “ 1-3x10'V3/4  L3/2  B2  (1-m2)  (251) 
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where  is  an  assumed  broad-hand  wave  amplitude.  For  D^  Schulz  and 
Lanzerotti  give  several  forms:  for  magnetic  impulses 
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and  for  electric  impulses 
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2.  SOLVING  THE  COMBINED  RADIAL  AND  PITCH-ANGLE  DIFFUSION  EQUATION 


The  combined  diffusion  equation  (248)  Is  well  suited  to  a numerical 
solution.  First  note  that  an  advantageous  change  of  variables  is 
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In  the  loss  cone  one  must  also  subtract  a loss  term 

/3f  \ 2wf 
UtL  Tb 

where  Tfa  is  the  bounce  period,  UR^T/V.  The  empirical  factor  w 
adjusted  to  match  the  detailed  solutions.  With  no  atmospheric 
scatter,  its  value  is  near  1 (Ref  41).  The  atmospheric  backscatter  at 
high  energies  is  slight  (Ref.  42);  it  is  probably  preferable  at  present  to 
ignore  backscatter  and  use  the  above  loss  rate  with  w=l.  The  form  of  the 
loss  term  is  such  that  one  need  only  replace  f by 


(259) 

can  be 
back- 


fLC  = 8 exp 


(260) 


and  solve  Eq.  (256)  (with  no  loss  term)  for  g in  place  of  f. 


Eq.  (256)  calls  for  a predictor-corrector  method  to  do  the  forward  time 
integration.  Let  f ' and  f'  be  the  x and  y diffusion  terms  on  the  right, 
thus 
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(261) 


An  Adams-Bashforth  type  method  is  suited  to  this  form.  The  essence  of 
that  method  is  to  expand  the  derivatives,  (~)  in  a Tayler  series. 


f'  (t)  » f'  (t  ) + 2C  , (t  - t ) + 3C  . (t  - t )2  + 


xl 


x2 


(262a) 


f’  (t)  * f'  (t  ) + 2C  . (t  - t ) + 3C  - (t  - t )2  + ...  (262b) 

y y n yl  n y2  n 

The  values  of  f and  t"  on  the  right  are  evaluated  at  the  point  t . 

x y 6 n 

The  predictor  step  consists  in  stepping  forward  from  the  point  t to  tn+^5 
using  the  values  of  f at  t^  and  the  values  of  f at  preceding  points 
(see  Figure  47).  The  predictor  step  is 
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In  the  corrector  step,  the  new  values  of  f at  cn+^  are  used  to  find 

f'  and  f'  there;  the  new  f^'s  are  then  used  to  find  new  C's  expanded 
x v 

about  t , The  corrector  step  then  gives  for  f: 
n+  L 
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r expansion  gives,  about  t : 

C”‘k-Ei  D.k 


(264) 


JL  , i vn-k 


mk  m+1 


(-1)  Klflall  A.^k,  ; 

.1  n 


(n-m-1)  at  a time  j > /TI  [all  A^,J;  m^n;k^n 


(265) 
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where ^ and  JT are  the  conventional  notation  for  extended  sums  and  products, 
respectively. 


The  formation  of  f'  and  f'  presents  difficulties  because  of  restrictions 
x y 

on  the  size  of  the  four-dimensional  integration  grid  (in  x,  y,  t,  and  n ) , 

It  pays  to  have  a differentiation  method  that  minimizes  computer  storage 
space,  even  at  the  apparent  expense  of  computing  time.  Spline  methods  are 
usually  quite  successful  in  estimating  derivatives  in  a grid  of  discrete 
points  (Ref.  43).  To  construct  a cubic  spline  of  f over  a discrete  variable 
£,  let 
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To  evaluate  the  derivatives,  (— ) ^ one  now  must  Invert  a tridiagonal 
matrix  equation 


Vi 


+ Vi+1 
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(271) 


The  spline,  however,  contains  two  unknown  parameters,  which  are  usually  to 
be  specified  at  the  two  end  points.  Derivative  boundary  conditions  are 
especially  useful,  because  two  (^)'s  are  directly  specified.  When  no 
boundary  condition  is  available,  it  is  necessary  to  estimate  the  derivative 
by  fitting  a polynomial  over  the  adjacent  end  points. 


Using  the  logarithm  in  place  of  f has  several  advantages,  the  most  important 
of  which  is  the  assurance  that  f will  remain  positive.  The  diffusion  opera- 
tors are  now  of  the  form 


(it  might  sometimes  be  desirable  to  compute  the  derivatives  of  f;  this 
is  rather  inefficient  because  f must  be  evaluated  from  log  (f)  at  every 
time  step.)  Two  conservative  applications  of  the  spline  yield  the  second 
order  differential  operators.  Figures  48  and  49  show  the  results  of 
several  trial  calculations.  A rather  coarse  grid  of  13  points  was  used,  none- 
theless the  second  derivatives  were  computed  quite  successfully.  A secondary 
advantage  of  Eq.  (272)  is  that  the  squared  term  on  the  right  usually 
dominates  for  interior  points.  This  aids  the  accuracy  in  minimizing  the 
errors  that  inevitably  occur  in  the  second  derivative. 


The  natural  boundary  conditions  for  the  radial  diffusion  operator  are  — **  0 
at  y 'v  0 and  f 0 at  f 'v  1.  The  first  condition  ensures  that  the  current 
(the  bracketed  term  in  Eq.  (256)  goes  to  zero  at  L ^ The  boundary  con- 
ditions for  the  pitch  angle  diffusion  operator  are  more  subtle.  If 

were  finite  at  the  ends,  x=0  and  x*l,  the  appropriate  boundary  condition 
3 f 

would  be  =0.  But  D may  vanish  at  the  ends,  leaving  an  indeterminate 

dX  XX 

boundary  condition.  In  that  case  one  obtains,  instead,  a relation  between 
the  derivatives: 
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Figure  48.  A test  case  showing  the  derivatives  evaluated 

by  a spline  approximation.  The  circled  points  are 
the  discrete  points  at  which  the  distribution  function, 
f,  was  evaluated.  This  case  corresponds  to  a high  L 
shell  with  a very  small  loss  cone  (on  the  right).  The 
exponents  denote  powers  of  10. 
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Figure  49.  A test  case  showing  the  derivatives  evaluated  by  a 

spline  approximation.  The  circled  points  are  the  discrete 

points  at  which  the  distribution  function,  f,  was  evaluated. 

8 f 

This  example  shows  clearly  the  excessive  curvature  of  n— 

o X 

near  x=0  which  makes  the  logarithmic  differentiation 
preferable.  This  case  is  for  an  intermediate  L shell.  The 
exponents  denote  powers  of  10. 


163 


» 

I 


l 


3_ 

3x 


il 

3x 


(273) 


3 f 

A supplementary  condition  may  also  be  posed  that  — is  always  negative 
at  the  two  ends;  this  ensures  that  numerical  errors  do  not  introduce  a 
spurious  particle  source  at  the  boundary. 

A computer  code  has  been  developed  to  implement  the  differential  equation 
solution.  It  constructs  solutions  over  constant  n slices  through  the 
distribution.  The  number  of  q’s  is  somewhat  arbitrary,  they  can  be  selected 
one  at  a time  and  the  results  can  be  stored  temporarily  so  that  the  resulting 
table  of  f values  can  be  interpolated  to  selected  energies.  The  main  deter- 
minants of  the  needed  storage  space  are  the  size  of  the  x,y  grid  and  the 
number  of  time  points  needed  for  a satisfactory  solution.  For  convenience 
the  x and  y arrays  are  nearly  the  same  size,  about  15  to  30  points  each. 

A satisfactory  integration  can  be  carried  out  with  less  than  5 to  8 time 
points.  The  maximum  number  of  f values  that  need  to  be  stored  for  any 
is  about  2000  (two  time  points).  The  maximum  number  of  f values  is  about 
15,000  (including  both  f'  and  f')»  with  approximately  the  same  number  of 
C's. 

In  the  cases  that  have  been  tested  there  were  always  regions  in  x,  y,  and  t 
where  radial  diffusion  dominated.  This  is  because,  despite  the  usual  vast 
differences  in  the  magnitude  of  the  diffusion  coefficients,  the  pitch  angle 
distribution  quickly  approaches  a lowest  normal  mode.  Steep  gradients  in  L, 
though,  can  persist  for  rather  long  times.  If  appreciable  differences  do 
exist  in  the  characteristic  diffusion  rates  the  scheme  developed  here  is 
quite  advantageous.  Usually  with  two  modes  of  diffusion  one  would  set  up 
two  separate  time  scales;  the  integration  method  described  above  is 
readily  adapted  to  that  approach.  We  selected  an  alternate  approach  that 
may  be  somewhat  less  efficient,  but  suffers  from  fewer  computational  complexi- 
ties. Our  approach  was  to  set  the  numbers  of  time  points  for  f^  and  f'  accord- 
ing to  the  respective  diffusion  rates.  The  use  of  a high  order  Adams- 
Bashforth  method  is  intended  to  maximize  the  time  steps.  The  step  size 
can  therefore  be  set  to  accommodate  the  slower  process,  with  just  2 or  3 
time  points;  the  order  of  the  faster  process  is  then  adjusted  for  acceptable 
accuracies.  Ideally  we  would  like  to  use  a vairable  order  at  each  point 
in  the  x,y  grid.  The  economies  of  such  a scheme  are  obvious,  although  cod- 
ing complexities  made  it  undesirable  for  the  first  trials. 
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3.  RESULTS  OF  RADIAL  DIFFUSION  COMPUTATIONS  AND  ADJUSTMENTS  TO  THE 
SPECTER  LOSS  MODEL 

In  the  earlier  computations  of  radial  diffusion  rates  following  an  artificial 
injection  event  (Ref.  35),  the  effects  of  radial  diffusion  were  negligible 
compared  with  loss  due  to  pitch-angle  diffusion.  Figure  50  is  a reproduc- 
tion of  Figure  3-3  from  the  earlier  report.  The  new  computations,  with 
increased  pitch  angle  diffusion  rates  at  early  times,  were  not  expected 
to  result  in  significantly  more  radial  diffusion  at  high  L-shells.  They 
should,  however,  tell  whether  the  elevated  radial  diffusion  rates  of  equa- 
torially  mirroring  electrons,  as  discussed  in  VI. 1,  could  alter  the  diffu- 
sion toward  low  L shells. 

Table  4 shows  the  results  of  the  first  few  steps  of  a trial  computation 
up  to  1.4  seconds. 

The  initial  distribution  was  chosen  to  have  a steep  radial  gradient  near 

L=3  and  a near-normal  mode  pitch-angle  distribution.  The  initial  values 

of  f'  and  f ' show  clearly  the  regions  which  favor  the  separate  diffusion 
x y 

processes.  The  example  does  demonstrate  the  expected  preference  for  radial 
diffusion  at  large  pitch  angles,  but  not  enough  to  outweigh  the  loss  over 
the  entire  distribution  due  to  pitch-angle  diffusion.  In  fact,  because  of 
the  different  radial  dependence  of  radial  and  pitch-angle  diffusion  coef- 
ficients, the  radial  gradient  in  the  omnidirectional  flux  can  actually 
steepen  with  time  at  the  outer  edge  of  the  "slot"  region  - L ^ 2.3  to 
L ffc  3.  A very  large  radial  diffusion  coefficient  was  chosen  for  the  test 
case,  in  a realistic  case  pitch-angle  diffusion  would  dominate  everywhere 
below  L % 6. 

It  may  be  concluded  that  radial  diffusion  is  not  likely  to  be  important  in 
artificial  radiation  belts,  except  perhaps  near  synchronous  altitudes.  At 
synchronous  altitudes  and  beyond  our  knowledge  of  radial  diffusion  rates 
is  not  sufficient  to  make  definitive  predictions.  Below  L •v#  3,  in  the 
slot  and  inner  belt,  radial  diffusion  can  be  definitely  ignored. 
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’igure  5o.  Jo>  "the  radial  omnidirectional  flux  at  the  equator,  at  various  times  for  the  detonation 
near  L = 4.0.  Shown  are  the  initial  flux  vaJ  for  t = 0 (solid  curve);  the  flux 
after  7 days  of  diffusion,  but  no  time  -e.-y.  determined  by  the  numerical  solution 
to  the  diffusion  equation  (dot-dashed  cur.  ):  d the  flux  after  7 days  of  time  decay, 

but  no  diffusion  (dashed  curve). 


Table  ^ Initial  results  of  a test  case  20.  The  unshaded  values  of  f^  and  f'  indicate  where 

each  dominates;  at  the  right  side  atmospheric  absorption  is  the  predominant  effect  Xn  particles  in  the 
loss  cone.  The  values  of  f(t=1.4)  are  shaded  where  a decrease  occurs.  The  diffusion  coefficients  at 
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The  calculations  of  radial  diffusion  rates  have  all  been  done  with  assumed 
rates  derived  for  the  natural  radiation  belts.  We  have  ignored  the  effect 
of  electromagnetic  disturbances  due  to  the  nuclear  explosion.  One  might, 
a priori,  conclude  that  the  disturbances  do  not  last  long  enough  to  cause 
substantial  radial  diffusion,  which  usually  takes  place  on  a time  scale  com' 
parable  to  the  particles  drift  periods.  The  discrepancies  noted  above 
between  the  injection  model  and  the  Telstar  data  above  L % 2 could  perhaps 
be  a result  of  radial  diffusion.  Whether  radial  diffusion  was  actually 
effective  is  a question  that  is  unlikely  to  be  answered.  Whether  the  dis- 
crepancies are  consistent  with  radial  diffusion  is  a question  that  can  be 
examined  with  the  new  diffusion  code. 


SECTION  VII 


RECONCILIATION  OF  THE  ELECTRON  LOSS  MODEL  TO  DATA  BODIES 


1.  NEW  DATA  BODIES 

In  an  earlier  report  (Reference  3)  an  extensive  comparison  was  made  of  the  Specter 
electron  loss  model  with  available  data,  particularly  from  the  artificial 
radiation  belts.  The  data  on  artificial  belts  have  been  quite  thoroughly  analyzed, 
and  no  high  altitude  nuclear  tests  are  expected  in  the  future,  so  we  must  rely  on 
natural  belt  data  to  fill  the  gaps  in  our  knowledge  of  decay  rates.  The  principal 
deficiencies  in  the  loss  model  are  at  low  altitudes,  where  the  longitude- 
dependent  filling  of  the  drift  loss  cone  is  important,  and  at  very  high  altitudes, 
where  the  decay  is  complicated  by  radial  diffusion  and  other  transient  phenomena. 

Several  important  data  bodies  have  been  reported  since  the  earlier  report 
(Reference  3).  At  low  and  intermediate  altitudes  the  relevant  new  data  are  those 
of  Imhof  (Reference  30)  and  Vampola  (Reference  29) . These  data  are  highly 
promising  since  they  often  define  the  loss  cone  distribution  to  a much  higher 
degree  of  precision  than  previous  data.  Unfortunately,  they  still  await  a detailed 
theoretical  analysis  that  might  define  the  rates  of  diffusion  into  the  loss  cone. 
This  problem  can  only  be  approached  with  an  extensive  computation  of  the  combined 
effects  of  pitch-angle  diffusion  at  high  altitudes,  scattering  in  the  atmosphere, 
and  longitudinal  drift.  Earlier  treatments,  though  marginally  adequate  for  the 
present  loss  model,  are  totally  insufficient  for  a precise  determination  of 
diffusion  rates,  mainly  because  of  neglect  of  atmospheric  backscatter.  Another 
problem  has  become  apparent  over  the  last  two  years,  to  wit:  what  is  the  signifi- 
cance of  interactions  with  man-made  VLF  emissions?  The  issue  is  still  somewhat 
controversial,  but  there  is  a growing  body  of  evidence  that  the  man-made  emissions 
are  very  important  (References  27  - 30).  At  high  altitudes  the  nature  of  the  VLF 
wave  sound  makes  little  difference  to  the  Specter  loss  model.  At  low  altitudes, 
because  of  the  complicating  effect  of  longitudinal  drift,  the  mechanism  of  wave 
generation  is  crucial.  Sporadic  or  coherent  VLF  waves  affect  trapped  electrons  in 


an  entirely  different  fashion  (Reference  31)  from  the  steady  loss  assumed  in  the 
theories  of  Kennel  and  Petscheck  (Reference  24)  and  Lyons,  Thorne,  and  Kennel 
(Reference  25).  The  issue  has  not  been  resolved  yet,  though  the  next  year  might 
bring  substantial  progress. 

At  high  altitudes,  the  recent  data  of  West  (Reference  44)  using  observations  from 
0G0-5,  are  especially  valuable.  For  the  first  time  we  have  available  data  with 
sufficient  resolution  to  test  the  combined  effects  of  radial  and  pitch  angle 
diffusion  beyond  L-4.  The  data  cover  the  energy  range  of  interest  in  Specter,  so 
their  analysis  should  yield  much  improved  loss  rates  in  the  outer  radiation  belt. 

The  data  are  just  beginning  to  be  analyzed  to  determine  diffusion  rates; 
unfortunately  no  concrete  results  were  available  at  the  time  of  this  report. 

2.  COMPARISON  OF  DATA  BODIES  WITH  SPECTER 

The  loss  rates  in  Specter  have  been  compared  with  all  the  available  artificial 
radiation  belt  data  and  large  parts  of  the  natural  radiation  belt  data,  the 
results  were  described  in  the  earlier  report  on  the  improvement  of  Specter  II 
(Reference  3).  The  only  major  unresolved  discrepancies  are  those  that  appear 
in  the  early  times,  as  discussed  above  in  Section  III.  Because  of  the  incomplete 
state  of  the  analysis  of  recent  data  we  were  not  able  to  make  any  additional  comparisons 
with  loss  rates  at  late  times. 

3.  IMPROVEMENTS  NEEDED  IN  THE  SPECTER  LOSS  MODEL;  NORMAL-MODE  DECAY  VS. 

NUMERICAL  MODELLING 

The  deficiencies  in  the  Specter  electron-loss  model  affect  the  decay  rate  at 
early  times  and  the  distribution  at  low  altitudes.  The  present  decay  model  is 
completely  empirical  and  obviously  cannot  apply  exactly  to  every  possible  case. 

A question  that  has  frequently  been  brought  up  is  the  desirability  of  replacing 
part  or  all  of  the  decay  model  by  an  actual  solution  to  a diffusion  equation. 

There  are  many  ways  this  might  be  accomplished,  among  them  the  use  of  an 
eigenfunction  expansion,  or  a direct  numerical  solution  of  a diffusion  equation, 
as  described  above  in  Section  VI. 
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The  use  of  an  eigenfunction  expansion  appears  attractive  because  of  the  possibility 
of  high  computational  efficiency.  Schulz  has  developed  an  interesting  scheme  in 
which  the  pitch-angle  variable  is  replaced  by  a new  variable  which  contains  the 
bounce  integral  and  all  the  sine  and  cosine  factors  taht  usually  appear  in  the 
diffusion  operator  (Reference  45).  His  variable  is  similar  to  x,  as  used  in 
Section  VI.  The  inclusion  of  the  bounce  integral,  T,  leads  to  inconsequential 
numerical  differences.  Schulz's  diffusion  coefficient  is,  however,  quite  different 
from  the  forms  assumed  above  for  Dxx.  The  diffusion  coefficient  DaC(  in  Section  VI 
was  assumed  finite,  which  led  to  awkward,  but  tractable,  boundary  conditions. 
Schulz,  on  the  other  hand,  assumed  a diffusion  coefficient  that  implied  Daa~°°  at 
both  ends,  x = 0°  and  a = 90°.  The  physical  implications  of  an  infinite  diffusion 
coefficient  are  not  clear;  the  published  values  (Reference  25)  of  Daa  seem 
relatively  constant  at  a~o  and  generally  vanish  at  a '90°.  The  Schulz  scheme  has 
one  special  advantage,  that  the  eigenfunctions  can  be  constructed  of  elemetnary 
functions  that  need  not  be  recomputed  at  every  step.  A loss  model  that  required 
the  evaluation  of  other  functions,  such  as  Bessel  functions,  would  be  out  of  the 
question  because  of  the  excessive  computation  time.  A criticism  that  may  be  very 
serious,  is  that  the  eigenfunctions  may  be  too  simple  to  adequately  represent  the 
lowest  normal  mode.  The  pitch-angle  distribution  of  the  trapped  particles  in 
Specter  was  carefully  fitted  to  a large  body  of  data  at  late  times.  It  is 
doubtful  whether  simpler  functions  would  be  adequate. 


An  alternate  approach  would  be  to  construct  a lowest  normal  mode  similar  to  the 
present  pitch-angle  distribution.  From  the  lowest  mode  the  higher  eigenfunctions 
could  be  computed  numerically  and  stored  as  interpolation  coefficients.  It  is 
quite  feasible  to  add  pitch-angle  points  (or  minor  points)  to  the  present  arrays 
to  facilitate  efficient,  accurate  computation.  Interpolation  by  evaluation  of  a 
second  or  third  degree  interpolation  polynomial  is  much  more  efficient  than  the 
evaluation  of  higher  functions.  (The  computation  of  the  eigenfunctions  is 
straightforward  and  need  only  be  performed  once,  so  the  added  computation  is 
Insignificant. ) 

The  (bounce  averaged)  pitch-angle  diffusion  equation  is 
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where  Q is  a source  of  particles.  In  a steady  state,  3f/8t=0,  a simple  solution 
results  for  D and  Q of  the  assumed  forms 
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The  above  diffusion  coefficient  is  physically  realistic  and  resembles  the 
coefficients  of  Lyons,  Thorne  and  Kennel  (Reference  25).  The  trapped 


distribution  is  then 


fa(v,2_y2+lni— + C ) 
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Where  uc  is  the  (global)  atmospheric  cutoff.  The  constant  C*  can  be  evaluated 
by  constructing  a loss  cone  solution  (involving  Bessel  functions)  (References 
46,  47)  and  fitting  the  solution  at  uc. 

Equation  (277)  suggests  a form  for  the  lowest  normal  mode  in  a freely-decaying 
distribution 


f = C, 


£(Wc“P2)n^^^ln  + Cj  J expf-t/T)  2<n 


(278) 


This  is  nearly  indistinguishable  from  the  distribution  now  in  Specter.  It  also 
has  the  advantage  of  a gradual  transition  across  yc  to  the  loss-cone  distribution. 
That  transition,  where  the  present  distribution  function  in  Specter  vanishes, 
has  been  a weak  point  in  the  representation  of  drift  loss-cone  fluxes.  Equation 
(278)  must  be  a solution  of  the  diffusion  equation 
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Substitution  in  (279)  gives  a differential  equation  for 
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Eq.  (280)  is  directly  integrable,  giving 
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Near  y ^ 0,  X must  approach  zero;  an  approximate  solution  is 
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This  gives 
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For  n=2 
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The  total  loss  rate  for  the  above  distribution  must  be  consistent  with  the 
rate  of  filling  of  the  loss  cone.  The  total  number  of  particles  contained 
ln  a tube  of  unit  cross  sectional  area  at  the  equator  is 


* - 2*  / f 4Rj.LT  ydy 


(285) 


'#  '<•! Iowa  f roe  (279)  and  (285)  that  the  total  loss  rate  is 
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f = 2*i  C ft  4relt^ 
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Substituting  for  D and  T gives 
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it  is  readily  seen  that  distributions  of  the  form 
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are  not  physically  realistic  approximations  to  a normal  decay  mode.  They 
should  be  replaced  by  the  more  complicated  formula,  (277),  to  ensure  that  the 
loss  cone  filling  is  computed  correctly.  Even  in  the  current  decay  model, 
this  is  an  important  consideration  because  the  precipitating  fluxes  at  low 
altitudes  are  normalized  to  rather  than  f. 

Replacing  the  present  long  time  distribution  by  the  above  distribution  (277) 
would  have  slight  effect  on  computational  efficiency,  and  would  ensure  that 
the  distribution  does  approach  a true  normal  mode.  This  distribution  could 
be  used  with  the  present  early  time  decay  model,  or,  alternately,  the  complete 
set  of  eigenfunctions  could  be  built  up  from  the  lowest  mode.  The  remaining 
eigenfunctions  could  be  easily  computed  with  the  aid  of  standard  eigenfunction 
generating  codes. 

Another,  quite  different,  approach  to  the  electron  decay  problem  is  to  employ 
a purely  numerical  method  that  works  for  any  arbitrary  set  of  diffusion 
coefficients.  An  appropriate  method  has  been  developed  and  described  above 
in  Section  VI.  The  computer  code  there  has  an  advantage  in  that  radial 
diffusion  is  explicitly  included;  though  radial  diffusion  can  be  omitted  with 
substantial  improvements  in  computing  efficiency.  The  code  at  present  is  not 
very  fast,  but  it  should  be  competitive  with  eigenfunction  methods.  Of 
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course,  the  explicit  integration  must  be  turned  off  when  it  reaches  the  lowest 
mode  decay  if  the  potential  computing  efficiency  is  to  be  realized.  All 
the  decay  models  have  one  feature  in  common:  that  they  revert  eventually  to 
a fixed  pitch-angle  distribution  which  can  be  either  represented  by  an 
algebraic  function,  or  by  a table  of  interpolation  coefficients. 

One  of  the  essential  features  of  any  decay  scheme  is  that  it  be  compatable 
with  the  model  used  to  derive  fluxes  in  the  drift  loss  cone.  It  has  not 
been  feasible  to  model  the  loss  cones  as  well  as  the  trapped  part  of  the 
distribution.  Perhaps  we  should  consider  a more  elaborate  calculation  in 
which  the  loss  cones  are  treated  on  an  equal  footing  to  the  rest  of  the 
distribution.  We  have  constructed  numerical  computations  of  the  drift  loss 
cone  filling  and  used  them  to  develop  the  current  loss  model.  That  calcula- 
tion did  not  treat  the  early-time  decay,  but  there  is  no  reason  why  all  phases 
of  the  decay  could  not  be  treated.  The  only  difficulty  is  in  modelling  the 
atmospheric  backscatter  (which  is  ignored  in  the  current  model).  We  now 
know  how  to  treat  the  atmospheric  backscatter,  and  approximate  it  quite  well 
in  a bounce  averaged  diffusion  calculation  (Refs.  11,  12).  A normal- 
mode  approach  may  be  feasible  for  the  drift  loss  cone  problem;  we  do  not, 
however,  know  of  any  way  the  atmospheric  backscatter  could  be  accurately 
represented  in  the  construction  of  eigenfunctions. 
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SECTION  VIII 

RECOMMENDATIONS  FOR  FUTURE  WORK 


t 


The  new  injection  model — including  the  initial  distribution  of  the  debris 
and  hot  air  near  the  magnetic  bubble,  the  initial  jetting  of  the  debris, 
the  outward  convection  of  the  debris  tube,  and  the  redistribution  of  the 
electrons  due  to  the  motion  of  the  debris  tube — is  a considerable  improve- 
ment over  the  previous  models.  The  pertinent  phenomena  are  treated  real- 
istically for  all  L values.  It  is  therefore  expected  to  provide  a good 
basis  for  the  computation  of  the  trapped-electron  distribution  due  to  any 
size  burst  at  any  point  in  the  magnetosphere,  not  just  for  the  limiting 
cases  of  the  very  small  and  very  large  bursts  at  points  where  the  tests 
have  been  conducted,  which  is  strictly  true  for  the  old  injection  models. 
Certain  parameters  of  the  old  models  had  been  adjusted  to  generate  fluxes 
that  were  in  agreement  with  the  nuclear  test  results.  There  are  no  such 
adjustable  parameters  in  the  new  model. 

Nevertheless,  the  model  still  has  shortcomings  that  should  be  improved 
in  future  efforts.  In  the  first  place,  the  assumption  in  the  debris  motion 
model  that  the  magnetic  field  lines  are  equipotentials  is  not  valid  for 
high  yields.  The  rate  of  change  of  the  neutralizing  current  along  the 
magnetic  field  is  so  large  that  the  electric  field  along  the  magnetic  field, 
due  to  the  effective  inductance  of  the  circuit,  is  not  negligible.  In 
fact,  for  Starfish  the  associated  potential  drop  along  the  magnetic  field 
at  early  times  was  about  one-third  of  the  potential  across  the  tube.  This 
effect  could  be  improved  by  using  a transmission-line  model,  with  approp- 
riate distributed  capacitance,  inductance,  and  resistance,  rather  than 
the  equivalent  LCR  resonant  circuit  with  lumped  parameters,  to  compute 
the  distribution  of  the  debris  in  the  magnetosphere.  In  addition  to  ac- 
counting for  the  potential  difference  along  the  magnetic  field,  this  approach 
would  also  correct  the  difficulty  in  the  present  model  that  the  entire 
capacitance  of  the  tube  is  used  in  the  circuit  even  while  the  hot-plasma 
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front  is  moving  along  the  tube.  This  effect  is  not  realistic  since  the 
thermal  plasma  ahead  of  the  front  does  not  know  about  the  advancing  plasma; 
i.e.,  the  front  moves  faster  than  an  Alfv4n  wave. 

As  discussed  in  Section  II. 2b,  subsection  (3),  it  appears  that  anomalous 
resistivity  may  also  appear  along  magnetic  field  lines.  Even  though  its 

t magnitude  cannot  be  determined  reliably  at  the  present  time,  represents- 

i 

| tive  values  based  on  observations  in  the  auroral  zones  could  be  included 

in  the  transmission-line  model  to  assess  its  importance  on  the  debris  dist- 
ribution. 

Futher  work  on  the  redistribution  model  is  also  warranted,  especially  since 
the  Telstar  satellite  data  obtained  after  the  Starfish  test  have  not  yet 
been  satisfactorily  explained  (see  Section  III. 2).  In  particular,  the  re- 
distribution of  the  electrons  that  occurs  during  their  eastward  drift  motion, 
while  they  are  bunched  in  longitude,  should  be  computed.  This  process  is 
discussed  in  Section  II. 4b. 

A problem  of  secondary  concern  is  the  likelihood  of  anomalous  trapping, 
as  originally  suggested  by  Davidson  (Ref.  48).  Though  the  expected  mag- 
nitude of  the  effect  is  not  very  large  (efficiencies  less  than  25%  have  been 
suggested  for  low  yield  bursts,  less  for  high  yield)  there  may  be  circum- 
stances, e.g.  nuclear  testing  at  80  to  200  kilometers  altitude,  where  there 
is  a possibility  of  serious  damage  to  satellites.  The  problem  merits  an 
exploratory  investigation  to  determine  its  likely  magnitude  and  impact  on 
SPECTER. 

A troublesome  shortcoming  of  the  trapped  electron  model  is  the  difficulty 
of  interpolating  meaningful  flux  values  to  obtain  the  energy  spectrum  at 
early  times.  As  discussed  in  Section  II. 5,  the  early-time  drift  diluted 
fluxes  are  highly  irregular  in  time,  longitude,  and  energy.  The  present 
interpolation  scheme  was  fully  adequate  as  long  as  the  total  flux  (or  1 
MeV  equivalent  flux)  alone  was  used  to  evaluate  satellite  irradiation. 

The  growing  need  for  spectral  information  points  out  the  desirability  of 
replacing  the  interpolation  of  drift  diluted  electron  fluxes.  The  most 
promising  approach  is  the  use  of  a semi-analytical  model,  suggested  in 
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Section  II. 5,  which  treats  the  drift  in  a fashion  consistent  with  first- 
principal  physics.  This  has  an  added  advantage  in  that  it  would  retain 
much  of  the  present  coding. 

Finally, the  computed  fluxes  for  Argus-3,  Starfish,  and  U.S.S.R.  tests 
should  be  further  tested  against  reliable  experimental  data. 

The  electron  loss  model  shows  some  obvious  needs  for  improvement,  mainly 
where  the  observational  data  have  been  inadequate  to  define  the  loss  rates. 
A continuing  search  for  better  data  is  called  for,  but  even  more  important 
is  a continuing  effort  to  incorporate  the  data  in  the  loss  model.  The 
loss  model  appears  adequate  for  most  of  the  trapped  electrons,  though 
much  remains  to  done  in  establishing  the  physical  bases  for  energy  and 
L-dependent  loss  rates.  A moderate  scale  program  should  be  directed 
toward  evaluating  the  improvements  that  could  be  expected  from  an  improved 
loss  model;  it  would  be  very  worthwhile  to  compare  the  advantages  of 
normal-mode  models  vs.  diffusion  models.  Only  after  such  a study  would 
it  be  wise  to  consider  replacing  the  present  trapped  electron  model  with 
something  more  complicated. 

The  low  altitude  part  of  the  loss  model,  especially  the  drift  loss  cone 
part,  is  not  so  satisfactory.  At  least  one  effect,  atmospheric  backscatter, 
has  been  totaly  neglected;  though  it  is  possible  to  achieve  a backscatter 
albedo  as  high  as  10  to  50%  at  1 to  3 MeV  for  a fission  decay  spectrum. 

The  backscatter  contribution  must  be  evaluated  in  a manner  consistent  with 
the  results  of  Davidson  and  Walt  (Ref.  41,  42).  The  relevant  observational 
data  have  not  yet  been  analyzed  to  determine  the  relative  contributions 
of  diffusion,  backscatter,  and  drift;  but,  nonetheless,  this  would  be 
an  appropriate  time  to  consider  improving  the  low  altitude  loss  model  in 
SPECTER.  A straightforward  approach  would  be  to  make  slight  improvements 
in  the  trapped  electron  loss  model,  as  discussed  in  Section  VII. 3,  and  do 
a new  set  of  calculations  for  fluxes  in  the  drift  loss  cone.  The  trapped 
and  loss  cone  distributions  could  be  joined  in  a way  similar  to  that  now 
used;  but  with  a more  detailed  pitch  angle  distribution. 
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The  SPECTER  development  program  la  presently  supported  by  several 
related  programs  at  Lockheed  Palo  Alto  Research  Laboratory.  We  have 
developed  a computer  code  to  calculate  pitch-angle  diffusion  into 
the  drift  loss  cones,  and  hope  to  apply  it  to  the  analysis  of  low 
altitude  energetic-election  data.  The  development  of  the  combined 
radial  and  pitch-angle  diffusion  code,  discussed  in  Section  VI 
is  continuing  independent  of  the  immediate  SPECTER  needs.  We  are  en- 
gaged in  a cooperative  program  with  Dr.  H.  West  of  Lawrence  Livermore 
Laboratories  to  use  the  code  to  analyze  his  outer  belt  electron  data 
(Ref.  44)  to  derive  loss  rates. 
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